This is gsl-ref.info, produced by makeinfo version 4.8 from gsl-ref.texi. INFO-DIR-SECTION Scientific software START-INFO-DIR-ENTRY * gsl-ref: (gsl-ref). GNU Scientific Library - Reference END-INFO-DIR-ENTRY Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 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.2 or any later version published by the Free Software Foundation; with the Invariant Sections being "GNU General Public License" and "Free Software Needs Free Documentation", the Front-Cover text being "A GNU Manual", and with the Back-Cover Text being (a) (see below). A copy of the license is included in the section entitled "GNU Free Documentation License". (a) The Back-Cover Text is: "You have freedom to copy and modify this GNU Manual, like GNU software."  File: gsl-ref.info, Node: Maximum and Minimum values, Next: Median and Percentiles, Prev: Weighted Samples, Up: Statistics 20.7 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 20.8 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 20.9 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 4.2984 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 20.10 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 `http://pdg.lbl.gov/'.  File: gsl-ref.info, Node: Histograms, Next: N-tuples, Prev: Statistics, Up: Top 21 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 21.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 21.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 21.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 21.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 21.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 21.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 21.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 21.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 21.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 21.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 21.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 21.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 21.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 21.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 21.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 21.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 21.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 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 21.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 21.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 21.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_histogram_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 21.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 21.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 22 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 22.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 22.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 22.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 22.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 22.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 22.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 22.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 22.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 22.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 23 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 23.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 23.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 / N) \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, double * XL, 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 23.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, double * XL, 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. The following variables can be accessed through the `gsl_monte_miser_state' struct, -- 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.  File: gsl-ref.info, Node: VEGAS, Next: Monte Carlo Examples, Prev: MISER, Up: Monte Carlo Integration 23.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 VEGAS algorithm is highly configurable. The following variables can be accessed through the `gsl_monte_vegas_state' struct, -- Variable: double result -- Variable: double sigma These parameters contain the raw value of the integral RESULT and its error SIGMA from the last iteration of the algorithm. -- Variable: double chisq This parameter gives the chi-squared per degree of freedom for the weighted estimate of the integral. The value of CHISQ should be close to 1. A value of CHISQ 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. -- 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.  File: gsl-ref.info, Node: Monte Carlo Examples, Next: Monte Carlo Integration References and Further Reading, Prev: VEGAS, Up: Monte Carlo Integration 23.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 = %.1g 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, s->chisq); } while (fabs (s->chisq - 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 0.6%. The estimated error `sigma' is consistent with the actual error, and the computed result differs from the true result by about one standard deviation, plain ================== result = 1.385867 sigma = 0.007938 exact = 1.393204 error = -0.007337 = 0.9 sigma The MISER algorithm reduces the error by a factor of two, and also correctly estimates the error, miser ================== result = 1.390656 sigma = 0.003743 exact = 1.393204 error = -0.002548 = 0.7 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.386925 sigma = 0.002651 exact = 1.393204 error = -0.006278 = 2 sigma converging... result = 1.392957 sigma = 0.000452 chisq/dof = 1.1 vegas final ================== result = 1.392957 sigma = 0.000452 exact = 1.393204 error = -0.000247 = 0.5 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 23.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 24 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 24.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 24.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: number_of_iterations temperature x x-(*x0_p) Ef(x) 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 24.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 24.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, $4}' | 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 24.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: -106.28 35.89 Los Alamos ###final_city_coord: -106.62 35.12 Albuquerque ###final_city_coord: -107.84 35.15 Grants ###final_city_coord: -107.87 37.29 Durango ###final_city_coord: -108.58 37.35 Cortez ###final_city_coord: -108.74 35.52 Gallup ###final_city_coord: -112.07 33.54 Phoenix ###final_city_coord: -106.76 32.34 Las Cruces ###final_city_coord: -96.77 32.79 Dallas ###final_city_coord: -103.2 34.41 Clovis ###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 24.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 25 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. These functions are declared in the header file `gsl_odeiv.h'. * Menu: * Defining the ODE System:: * Stepping Functions:: * Adaptive Step-size Control:: * Evolution:: * ODE Example programs:: * ODE References and Further Reading::  File: gsl-ref.info, Node: Defining the ODE System, Next: Stepping Functions, Up: Ordinary Differential Equations 25.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_odeiv_system' datatype. -- Data Type: gsl_odeiv_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. `int (* jacobian) (double t, const double y[], double * dfdy, double dfdt[], void * params);' This function should store the vector of derivative elements df_i(t,y,params)/dt 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. The function should return `GSL_SUCCESS' if the calculation was completed successfully. Any other return value indicates an error. Some of the simpler solver algorithms do not make use of the Jacobian matrix, so it is not always strictly necessary to provide it (the `jacobian' element of the struct can be replaced by a null pointer for those algorithms). However, it is useful to provide the Jacobian to allow the solver algorithms to be interchanged--the best algorithms make use of the Jacobian. `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 25.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_odeiv_step * gsl_odeiv_step_alloc (const gsl_odeiv_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. -- Function: int gsl_odeiv_step_reset (gsl_odeiv_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_odeiv_step_free (gsl_odeiv_step * S) This function frees all the memory associated with the stepping function S. -- Function: const char * gsl_odeiv_step_name (const gsl_odeiv_step * S) This function returns a pointer to the name of the stepping function. For example, printf ("step method is '%s'\n", gsl_odeiv_step_name (s)); would print something like `step method is 'rk4''. -- Function: unsigned int gsl_odeiv_step_order (const gsl_odeiv_step * S) This function returns the order of the stepping function on the previous step. This order can vary if the stepping function itself is adaptive. -- Function: int gsl_odeiv_step_apply (gsl_odeiv_step * S, double T, double H, double Y[], double YERR[], const double DYDT_IN[], double DYDT_OUT[], const gsl_odeiv_system * DYDT) This function applies the stepping function S to the system of equations defined by DYDT, 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. If the user-supplied functions defined in the system DYDT return a status other than `GSL_SUCCESS' the step will be aborted. In this 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. The step-size H will be set to the step-size which caused the error. If the function is called again with a smaller step-size, e.g. H/10, it should be possible to get closer to any singularity. To distinguish between error codes from the user-supplied functions and those from `gsl_odeiv_step_apply' itself, any user-defined return values should be distinct from the standard GSL error codes. The following algorithms are available, -- Step Type: gsl_odeiv_step_rk2 Embedded Runge-Kutta (2, 3) method. -- Step Type: gsl_odeiv_step_rk4 4th order (classical) Runge-Kutta. -- Step Type: gsl_odeiv_step_rkf45 Embedded Runge-Kutta-Fehlberg (4, 5) method. This method is a good general-purpose integrator. -- Step Type: gsl_odeiv_step_rkck Embedded Runge-Kutta Cash-Karp (4, 5) method. -- Step Type: gsl_odeiv_step_rk8pd Embedded Runge-Kutta Prince-Dormand (8,9) method. -- Step Type: gsl_odeiv_step_rk2imp Implicit 2nd order Runge-Kutta at Gaussian points. -- Step Type: gsl_odeiv_step_rk4imp Implicit 4th order Runge-Kutta at Gaussian points. -- Step Type: gsl_odeiv_step_bsimp Implicit Bulirsch-Stoer method of Bader and Deuflhard. This algorithm requires the Jacobian. -- Step Type: gsl_odeiv_step_gear1 M=1 implicit Gear method. -- Step Type: gsl_odeiv_step_gear2 M=2 implicit Gear method.  File: gsl-ref.info, Node: Adaptive Step-size Control, Next: Evolution, Prev: Stepping Functions, Up: Ordinary Differential Equations 25.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_odeiv_control * gsl_odeiv_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'_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_odeiv_control * gsl_odeiv_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_odeiv_control * gsl_odeiv_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_odeiv_control * gsl_odeiv_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_odeiv_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'_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_odeiv_control * gsl_odeiv_control_alloc (const gsl_odeiv_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_odeiv_control_init (gsl_odeiv_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_odeiv_control_free (gsl_odeiv_control * C) This function frees all the memory associated with the control function C. -- Function: int gsl_odeiv_control_hadjust (gsl_odeiv_control * C, gsl_odeiv_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_odeiv_control_name (const gsl_odeiv_control * C) This function returns a pointer to the name of the control function. For example, printf ("control method is '%s'\n", gsl_odeiv_control_name (c)); would print something like `control method is 'standard''  File: gsl-ref.info, Node: Evolution, Next: ODE Example programs, Prev: Adaptive Step-size Control, Up: Ordinary Differential Equations 25.4 Evolution ============== The highest level of the system is the evolution function which combines the results of a stepping function and control function to reliably advance the solution forward over an interval (t_0, t_1). If the control function signals that the step-size should be decreased the evolution function backs out of the current step and tries the proposed smaller step-size. This process is continued until an acceptable step-size is found. -- Function: gsl_odeiv_evolve * gsl_odeiv_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_odeiv_evolve_apply (gsl_odeiv_evolve * E, gsl_odeiv_control * CON, gsl_odeiv_step * STEP, const gsl_odeiv_system * DYDT, double * T, double T1, double * H, double Y[]) This function advances the system (E, DYDT) 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, but this will be modified using the control function C to achieve the appropriate error bound if necessary. The routine may make several calls to STEP in order to determine the optimum step-size. An estimate of the local error for the step can be obtained from the components of the array `E->yerr[]'. If the step-size has been changed the value of H will be modified on output. 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. If the user-supplied functions defined in the system DYDT return a status other than `GSL_SUCCESS' the step will be aborted. In this case, T and Y will be restored to their pre-step values and the error code from the user-supplied function will be returned. To distinguish between error codes from the user-supplied functions and those from `gsl_odeiv_evolve_apply' itself, any user-defined return values should be distinct from the standard GSL error codes. -- Function: int gsl_odeiv_evolve_reset (gsl_odeiv_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_odeiv_evolve_free (gsl_odeiv_evolve * E) This function frees all the memory associated with the evolution function E.  File: gsl-ref.info, Node: ODE Example programs, Next: ODE References and Further Reading, Prev: Evolution, Up: Ordinary Differential Equations 25.5 Examples ============= The following program solves the second-order nonlinear Van der Pol oscillator equation, x''(t) + \mu x'(t) (x(t)^2 - 1) + x(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, y = x'(t), x' = y y' = -x + \mu y (1-x^2) The program begins by defining functions for these derivatives and their Jacobian, #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) { const gsl_odeiv_step_type * T = gsl_odeiv_step_rk8pd; gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, 2); gsl_odeiv_control * c = gsl_odeiv_control_y_new (1e-6, 0.0); gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (2); double mu = 10; gsl_odeiv_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_odeiv_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_odeiv_evolve_free (e); gsl_odeiv_control_free (c); gsl_odeiv_step_free (s); return 0; } For functions with multiple parameters, the appropriate information can be passed in through the PARAMS argument using a pointer to a struct. The main loop of the program evolves the solution from (y, y') = (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 Y. To obtain the values at regular intervals, rather than the variable spacings chosen by the control function, the main loop can be modified to advance the solution from one point to the next. For example, the following main loop prints the solution at the fixed points t = 0, 1, 2, \dots, 100, for (i = 1; i <= 100; i++) { double ti = i * t1 / 100.0; while (t < ti) { gsl_odeiv_evolve_apply (e, c, s, &sys, &t, ti, &h, y); } printf ("%.5e %.5e %.5e\n", t, y[0], y[1]); } It is also possible to work with a non-adaptive integrator, using only the stepping function itself. The following program uses the `rk4' fourth-order Runge-Kutta stepping function with a fixed stepsize of 0.01, int main (void) { const gsl_odeiv_step_type * T = gsl_odeiv_step_rk4; gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, 2); double mu = 10; gsl_odeiv_system sys = {func, jac, 2, &mu}; double t = 0.0, t1 = 100.0; double h = 1e-2; double y[2] = { 1.0, 0.0 }, y_err[2]; double dydt_in[2], dydt_out[2]; /* initialise dydt_in from system parameters */ GSL_ODEIV_FN_EVAL(&sys, t, y, dydt_in); while (t < t1) { int status = gsl_odeiv_step_apply (s, t, h, y, y_err, dydt_in, dydt_out, &sys); if (status != GSL_SUCCESS) break; dydt_in[0] = dydt_out[0]; dydt_in[1] = dydt_out[1]; t += h; printf ("%.5e %.5e %.5e\n", t, y[0], y[1]); } gsl_odeiv_step_free (s); return 0; } The derivatives must be initialized for the starting point t=0 before the first step is taken. Subsequent steps use the output derivatives DYDT_OUT as inputs to the next step by copying their values into DYDT_IN.  File: gsl-ref.info, Node: ODE References and Further Reading, Prev: ODE Example programs, Up: Ordinary Differential Equations 25.6 References and Further Reading =================================== 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.  File: gsl-ref.info, Node: Interpolation, Next: Numerical Differentiation, Prev: Ordinary Differential Equations, Up: Top 26 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. 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 26.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 26.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; 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 26.3 Interpolation Types ======================== The interpolation library provides five 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) This function returns the minimum number of points required by the interpolation type of INTERP. 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 26.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]. -- 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]'. -- 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 26.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. -- 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 26.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 26.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 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: first = last 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 26.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 27 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 27.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 27.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 27.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 28 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:: * 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 28.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: Chebyshev Series Evaluation, Prev: Chebyshev Definitions, Up: Chebyshev Approximations 28.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: Chebyshev Series Evaluation, Next: Derivatives and Integrals, Prev: Creation and Calculation of Chebyshev Series, Up: Chebyshev Approximations 28.3 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 N, 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 28.4 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 28.5 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 28.6 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 29 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 29.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 29.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 29.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 29.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, `http://arxiv.org/abs/math/0005209'.  File: gsl-ref.info, Node: Wavelet Transforms, Next: Discrete Hankel Transforms, Prev: Series Acceleration, Up: Top 30 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 30.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 30.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 The 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. -- 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 30.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 30.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 30.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 and non-standard forms 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 30.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 30.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 `http://www.physionet.org/' 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 31 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 31.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. 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)x) / 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 31.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 X. -- 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 X. -- 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. -- 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 31.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 32 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 32.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 32.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 32.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 32.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: 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 following macro, #define GSL_FN_EVAL(F,x) (*((F)->function))(x,(F)->params) -- 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 following macro, #define GSL_FN_FDF_EVAL_F(FDF,x) (*((FDF)->f))(x,(FDF)->params) The derivative f'(x) can be evaluated using the following macro, #define GSL_FN_FDF_EVAL_DF(FDF,x) (*((FDF)->df))(x,(FDF)->params) and both the function y = f(x) and its derivative dy = f'(x) can be evaluated at the same time using the following macro, #define GSL_FN_FDF_EVAL_F_DF(FDF,x,y,dy) (*((FDF)->fdf))(x,(FDF)->params,(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 32.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 a 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 32.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 32.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 32.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 32.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" 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.  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 32.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 32.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 33 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 33.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 33.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 33.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 33.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 33.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 33.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 33.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.  File: gsl-ref.info, Node: Minimization Examples, Next: Minimization References and Further Reading, Prev: Minimization Algorithms, Up: One dimensional Minimization 33.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 %.7f\n", iter, a, b, m, m_expected, 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  File: gsl-ref.info, Node: Minimization References and Further Reading, Prev: Minimization Examples, Up: One dimensional Minimization 33.9 References and Further Reading =================================== Further information on Brent's algorithm is available in the following book, Richard Brent, `Algorithms for minimization without derivatives', Prentice-Hall (1973), republished by Dover in paperback (2002), ISBN 0-486-41998-3.  File: gsl-ref.info, Node: Multidimensional Root-Finding, Next: Multidimensional Minimization, Prev: One dimensional Minimization, Up: Top 34 Multidimensional Root-Finding ******************************** This chapter describes functions for multidimensional root-finding (solving nonlinear systems with n equations in n unknowns). 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 solvers are based on the original Fortran library MINPACK. The header file `gsl_multiroots.h' contains prototypes for the multidimensional root finding functions and related declarations. * Menu: * Overview of Multidimensional Root Finding:: * Initializing the Multidimensional Solver:: * Providing the multidimensional system of equations to solve:: * Iteration of the multidimensional solver:: * Search Stopping Parameters for the multidimensional solver:: * Algorithms using Derivatives:: * Algorithms without Derivatives:: * Example programs for Multidimensional Root finding:: * References and Further Reading for Multidimensional Root Finding::  File: gsl-ref.info, Node: Overview of Multidimensional Root Finding, Next: Initializing the Multidimensional Solver, Up: Multidimensional Root-Finding 34.1 Overview ============= The problem of multidimensional root finding requires the simultaneous solution of n equations, f_i, in n variables, x_i, f_i (x_1, ..., x_n) = 0 for i = 1 ... n. In general there are no bracketing methods available for n dimensional systems, and no way of knowing whether any solutions exist. All algorithms proceed from an initial guess using a variant of the Newton iteration, x -> x' = x - J^{-1} f(x) where x, f are vector quantities and J is the Jacobian matrix J_{ij} = d f_i / d x_j. Additional strategies can be used to enlarge the region of convergence. These include requiring a decrease in the norm |f| on each step proposed by Newton's method, or taking steepest-descent steps in the direction of the negative gradient of |f|. Several root-finding algorithms are available within a single framework. 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 evaluation of the Jacobian matrix can be problematic, either because programming the derivatives is intractable or because computation of the n^2 terms of the matrix becomes too expensive. For these reasons the algorithms provided by the library are divided into two classes according to whether the derivatives are available or not. The state for solvers with an analytic Jacobian matrix is held in a `gsl_multiroot_fdfsolver' struct. The updating procedure requires both the function and its derivatives to be supplied by the user. The state for solvers which do not use an analytic Jacobian matrix is held in a `gsl_multiroot_fsolver' struct. The updating procedure uses only function evaluations (not derivatives). The algorithms estimate the matrix J or J^{-1} by approximate methods.  File: gsl-ref.info, Node: Initializing the Multidimensional Solver, Next: Providing the multidimensional system of equations to solve, Prev: Overview of Multidimensional Root Finding, Up: Multidimensional Root-Finding 34.2 Initializing the Solver ============================ The following functions initialize a multidimensional solver, either with or without derivatives. The solver itself depends only on the dimension of the problem and the algorithm and can be reused for different problems. -- Function: gsl_multiroot_fsolver * gsl_multiroot_fsolver_alloc (const gsl_multiroot_fsolver_type * T, size_t N) This function returns a pointer to a newly allocated instance of a solver of type T for a system of N dimensions. For example, the following code creates an instance of a hybrid solver, to solve a 3-dimensional system of equations. const gsl_multiroot_fsolver_type * T = gsl_multiroot_fsolver_hybrid; gsl_multiroot_fsolver * s = gsl_multiroot_fsolver_alloc (T, 3); 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_multiroot_fdfsolver * gsl_multiroot_fdfsolver_alloc (const gsl_multiroot_fdfsolver_type * T, size_t N) This function returns a pointer to a newly allocated instance of a derivative solver of type T for a system of N dimensions. For example, the following code creates an instance of a Newton-Raphson solver, for a 2-dimensional system of equations. const gsl_multiroot_fdfsolver_type * T = gsl_multiroot_fdfsolver_newton; gsl_multiroot_fdfsolver * s = gsl_multiroot_fdfsolver_alloc (T, 2); 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_multiroot_fsolver_set (gsl_multiroot_fsolver * S, gsl_multiroot_function * F, gsl_vector * X) This function sets, or resets, an existing solver S to use the function F and the initial guess X. -- Function: int gsl_multiroot_fdfsolver_set (gsl_multiroot_fdfsolver * S, gsl_multiroot_function_fdf * FDF, gsl_vector * X) This function sets, or resets, an existing solver S to use the function and derivative FDF and the initial guess X. -- Function: void gsl_multiroot_fsolver_free (gsl_multiroot_fsolver * S) -- Function: void gsl_multiroot_fdfsolver_free (gsl_multiroot_fdfsolver * S) These functions free all the memory associated with the solver S. -- Function: const char * gsl_multiroot_fsolver_name (const gsl_multiroot_fsolver * S) -- Function: const char * gsl_multiroot_fdfsolver_name (const gsl_multiroot_fdfsolver * S) These functions return a pointer to the name of the solver. For example, printf ("s is a '%s' solver\n", gsl_multiroot_fdfsolver_name (s)); would print something like `s is a 'newton' solver'.  File: gsl-ref.info, Node: Providing the multidimensional system of equations to solve, Next: Iteration of the multidimensional solver, Prev: Initializing the Multidimensional Solver, Up: Multidimensional Root-Finding 34.3 Providing the function to solve ==================================== You must provide n functions of n variables for the root finders to operate on. In order to allow for general parameters the functions are defined by the following data types: -- Data Type: gsl_multiroot_function This data type defines a general system of functions with parameters. `int (* f) (const gsl_vector * X, void * PARAMS, gsl_vector * F)' this function should store the vector result f(x,params) in F for argument X and parameters PARAMS, returning an appropriate error code if the function cannot be computed. `size_t n' the dimension of the system, i.e. the number of components of the vectors X and F. `void * params' a pointer to the parameters of the function. Here is an example using Powell's test function, f_1(x) = A x_0 x_1 - 1, f_2(x) = exp(-x_0) + exp(-x_1) - (1 + 1/A) with A = 10^4. The following code defines a `gsl_multiroot_function' system `F' which you could pass to a solver: struct powell_params { double A; }; int powell (gsl_vector * x, void * p, gsl_vector * f) { struct powell_params * params = *(struct powell_params *)p; const double A = (params->A); const double x0 = gsl_vector_get(x,0); const double x1 = gsl_vector_get(x,1); gsl_vector_set (f, 0, A * x0 * x1 - 1); gsl_vector_set (f, 1, (exp(-x0) + exp(-x1) - (1.0 + 1.0/A))); return GSL_SUCCESS } gsl_multiroot_function F; struct powell_params params = { 10000.0 }; F.f = &powell; F.n = 2; F.params = ¶ms; -- Data Type: gsl_multiroot_function_fdf This data type defines a general system of functions with parameters and the corresponding Jacobian matrix of derivatives, `int (* f) (const gsl_vector * X, void * PARAMS, gsl_vector * F)' this function should store the vector result f(x,params) in F for argument X and parameters PARAMS, returning an appropriate error code if the function cannot be computed. `int (* df) (const gsl_vector * X, void * PARAMS, gsl_matrix * J)' this function should store the N-by-N matrix result J_ij = d f_i(x,params) / d x_j in J for argument X and parameters PARAMS, returning an appropriate error code if the function cannot be computed. `int (* fdf) (const gsl_vector * X, void * PARAMS, gsl_vector * F, gsl_matrix * J)' This function should set the values of the F and J as above, for arguments X and parameters PARAMS. This function provides an optimization of the separate functions for f(x) and J(x)--it is always faster to compute the function and its derivative at the same time. `size_t n' the dimension of the system, i.e. the number of components of the vectors X and F. `void * params' a pointer to the parameters of the function. The example of Powell's test function defined above can be extended to include analytic derivatives using the following code, int powell_df (gsl_vector * x, void * p, gsl_matrix * J) { struct powell_params * params = *(struct powell_params *)p; const double A = (params->A); const double x0 = gsl_vector_get(x,0); const double x1 = gsl_vector_get(x,1); gsl_matrix_set (J, 0, 0, A * x1); gsl_matrix_set (J, 0, 1, A * x0); gsl_matrix_set (J, 1, 0, -exp(-x0)); gsl_matrix_set (J, 1, 1, -exp(-x1)); return GSL_SUCCESS } int powell_fdf (gsl_vector * x, void * p, gsl_matrix * f, gsl_matrix * J) { struct powell_params * params = *(struct powell_params *)p; const double A = (params->A); const double x0 = gsl_vector_get(x,0); const double x1 = gsl_vector_get(x,1); const double u0 = exp(-x0); const double u1 = exp(-x1); gsl_vector_set (f, 0, A * x0 * x1 - 1); gsl_vector_set (f, 1, u0 + u1 - (1 + 1/A)); gsl_matrix_set (J, 0, 0, A * x1); gsl_matrix_set (J, 0, 1, A * x0); gsl_matrix_set (J, 1, 0, -u0); gsl_matrix_set (J, 1, 1, -u1); return GSL_SUCCESS } gsl_multiroot_function_fdf FDF; FDF.f = &powell_f; FDF.df = &powell_df; FDF.fdf = &powell_fdf; FDF.n = 2; FDF.params = 0; Note that the function `powell_fdf' is able to reuse existing terms from the function when calculating the Jacobian, thus saving time.  File: gsl-ref.info, Node: Iteration of the multidimensional solver, Next: Search Stopping Parameters for the multidimensional solver, Prev: Providing the multidimensional system of equations to solve, Up: Multidimensional Root-Finding 34.4 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_multiroot_fsolver_iterate (gsl_multiroot_fsolver * S) -- Function: int gsl_multiroot_fdfsolver_iterate (gsl_multiroot_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_ENOPROG' the iteration is not making any progress, preventing the algorithm from continuing. The solver maintains a current best estimate of the root at all times. This information can be accessed with the following auxiliary functions, -- Function: gsl_vector * gsl_multiroot_fsolver_root (const gsl_multiroot_fsolver * S) -- Function: gsl_vector * gsl_multiroot_fdfsolver_root (const gsl_multiroot_fdfsolver * S) These functions return the current estimate of the root for the solver S. -- Function: gsl_vector * gsl_multiroot_fsolver_f (const gsl_multiroot_fsolver * S) -- Function: gsl_vector * gsl_multiroot_fdfsolver_f (const gsl_multiroot_fdfsolver * S) These functions return the function value f(x) at the current estimate of the root for the solver S. -- Function: gsl_vector * gsl_multiroot_fsolver_dx (const gsl_multiroot_fsolver * S) -- Function: gsl_vector * gsl_multiroot_fdfsolver_dx (const gsl_multiroot_fdfsolver * S) These functions return the last step dx taken by the solver S.  File: gsl-ref.info, Node: Search Stopping Parameters for the multidimensional solver, Next: Algorithms using Derivatives, Prev: Iteration of the multidimensional solver, Up: Multidimensional Root-Finding 34.5 Search Stopping Parameters =============================== A root finding procedure should stop when one of the following conditions is true: * A multidimensional 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_multiroot_test_delta (const gsl_vector * DX, const gsl_vector * X, double EPSABS, double EPSREL) This function tests for the convergence of the sequence by comparing the last step DX with the absolute error EPSABS and relative error EPSREL to the current position X. The test returns `GSL_SUCCESS' if the following condition is achieved, |dx_i| < epsabs + epsrel |x_i| for each component of X and returns `GSL_CONTINUE' otherwise. -- Function: int gsl_multiroot_test_residual (const gsl_vector * 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, \sum_i |f_i| < 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 is small enough.  File: gsl-ref.info, Node: Algorithms using Derivatives, Next: Algorithms without Derivatives, Prev: Search Stopping Parameters for the multidimensional solver, Up: Multidimensional Root-Finding 34.6 Algorithms using Derivatives ================================= The root finding algorithms described in this section make use of both the function and its derivative. They require an initial guess for the location of the root, but 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 the conditions are satisfied then convergence is quadratic. -- Derivative Solver: gsl_multiroot_fdfsolver_hybridsj This is a modified version of Powell's Hybrid method as implemented in the HYBRJ algorithm in MINPACK. Minpack was written by Jorge J. More', Burton S. Garbow and Kenneth E. Hillstrom. The Hybrid algorithm retains the fast convergence of Newton's method but will also reduce the residual when Newton's method is unreliable. The algorithm uses a generalized trust region to keep each step under control. In order to be accepted a proposed new position x' must satisfy the condition |D (x' - x)| < \delta, where D is a diagonal scaling matrix and \delta is the size of the trust region. The components of D are computed internally, using the column norms of the Jacobian to estimate the sensitivity of the residual to each component of x. This improves the behavior of the algorithm for badly scaled functions. On each iteration the algorithm first determines the standard Newton step by solving the system J dx = - f. If this step falls inside the trust region it is used as a trial step in the next stage. If not, the algorithm uses the linear combination of the Newton and gradient directions which is predicted to minimize the norm of the function while staying inside the trust region, dx = - \alpha J^{-1} f(x) - \beta \nabla |f(x)|^2. This combination of Newton and gradient directions is referred to as a "dogleg step". The proposed step is now tested by evaluating the function at the resulting point, x'. If the step reduces the norm of the function sufficiently then it is accepted and size of the trust region is increased. If the proposed step fails to improve the solution then the size of the trust region is decreased and another trial step is computed. The speed of the algorithm is increased by computing the changes to the Jacobian approximately, using a rank-1 update. If two successive attempts fail to reduce the residual then the full Jacobian is recomputed. The algorithm also monitors the progress of the solution and returns an error if several steps fail to make any improvement, `GSL_ENOPROG' the iteration is not making any progress, preventing the algorithm from continuing. `GSL_ENOPROGJ' re-evaluations of the Jacobian indicate that the iteration is not making any progress, preventing the algorithm from continuing. -- Derivative Solver: gsl_multiroot_fdfsolver_hybridj This algorithm is an unscaled version of `hybridsj'. The steps are controlled by a spherical trust region |x' - x| < \delta, instead of a generalized region. This can be useful if the generalized region estimated by `hybridsj' is inappropriate. -- Derivative Solver: gsl_multiroot_fdfsolver_newton Newton's Method is the standard root-polishing algorithm. The algorithm begins with an initial guess for the location of the solution. On each iteration a linear approximation to the function F is used to estimate the step which will zero all the components of the residual. The iteration is defined by the following sequence, x -> x' = x - J^{-1} f(x) where the Jacobian matrix J is computed from the derivative functions provided by F. The step dx is obtained by solving the linear system, J dx = - f(x) using LU decomposition. -- Derivative Solver: gsl_multiroot_fdfsolver_gnewton This is a modified version of Newton's method which attempts to improve global convergence by requiring every step to reduce the Euclidean norm of the residual, |f(x)|. If the Newton step leads to an increase in the norm then a reduced step of relative size, t = (\sqrt(1 + 6 r) - 1) / (3 r) is proposed, with r being the ratio of norms |f(x')|^2/|f(x)|^2. This procedure is repeated until a suitable step size is found.  File: gsl-ref.info, Node: Algorithms without Derivatives, Next: Example programs for Multidimensional Root finding, Prev: Algorithms using Derivatives, Up: Multidimensional Root-Finding 34.7 Algorithms without Derivatives =================================== The algorithms described in this section do not require any derivative information to be supplied by the user. Any derivatives needed are approximated by finite differences. Note that if the finite-differencing step size chosen by these routines is inappropriate, an explicit user-supplied numerical derivative can always be used with the algorithms described in the previous section. -- Solver: gsl_multiroot_fsolver_hybrids This is a version of the Hybrid algorithm which replaces calls to the Jacobian function by its finite difference approximation. The finite difference approximation is computed using `gsl_multiroots_fdjac' with a relative step size of `GSL_SQRT_DBL_EPSILON'. Note that this step size will not be suitable for all problems. -- Solver: gsl_multiroot_fsolver_hybrid This is a finite difference version of the Hybrid algorithm without internal scaling. -- Solver: gsl_multiroot_fsolver_dnewton The "discrete Newton algorithm" is the simplest method of solving a multidimensional system. It uses the Newton iteration x -> x - J^{-1} f(x) where the Jacobian matrix J is approximated by taking finite differences of the function F. The approximation scheme used by this implementation is, J_{ij} = (f_i(x + \delta_j) - f_i(x)) / \delta_j where \delta_j is a step of size \sqrt\epsilon |x_j| with \epsilon being the machine precision (\epsilon \approx 2.22 \times 10^-16). The order of convergence of Newton's algorithm is quadratic, but the finite differences require n^2 function evaluations on each iteration. The algorithm may become unstable if the finite differences are not a good approximation to the true derivatives. -- Solver: gsl_multiroot_fsolver_broyden The "Broyden algorithm" is a version of the discrete Newton algorithm which attempts to avoids the expensive update of the Jacobian matrix on each iteration. The changes to the Jacobian are also approximated, using a rank-1 update, J^{-1} \to J^{-1} - (J^{-1} df - dx) dx^T J^{-1} / dx^T J^{-1} df where the vectors dx and df are the changes in x and f. On the first iteration the inverse Jacobian is estimated using finite differences, as in the discrete Newton algorithm. This approximation gives a fast update but is unreliable if the changes are not small, and the estimate of the inverse Jacobian becomes worse as time passes. The algorithm has a tendency to become unstable unless it starts close to the root. The Jacobian is refreshed if this instability is detected (consult the source for details). This algorithm is included only for demonstration purposes, and is not recommended for serious use.  File: gsl-ref.info, Node: Example programs for Multidimensional Root finding, Next: References and Further Reading for Multidimensional Root Finding, Prev: Algorithms without Derivatives, Up: Multidimensional Root-Finding 34.8 Examples ============= The multidimensional solvers are used in a similar way to the one-dimensional root finding algorithms. This first example demonstrates the `hybrids' scaled-hybrid algorithm, which does not require derivatives. The program solves the Rosenbrock system of equations, f_1 (x, y) = a (1 - x) f_2 (x, y) = b (y - x^2) with a = 1, b = 10. The solution of this system lies at (x,y) = (1,1) in a narrow valley. The first stage of the program is to define the system of equations, #include #include #include #include struct rparams { double a; double b; }; int rosenbrock_f (const gsl_vector * x, void *params, gsl_vector * f) { double a = ((struct rparams *) params)->a; double b = ((struct rparams *) params)->b; const double x0 = gsl_vector_get (x, 0); const double x1 = gsl_vector_get (x, 1); const double y0 = a * (1 - x0); const double y1 = b * (x1 - x0 * x0); gsl_vector_set (f, 0, y0); gsl_vector_set (f, 1, y1); return GSL_SUCCESS; } The main program begins by creating the function object `f', with the arguments `(x,y)' and parameters `(a,b)'. The solver `s' is initialized to use this function, with the `hybrids' method. int main (void) { const gsl_multiroot_fsolver_type *T; gsl_multiroot_fsolver *s; int status; size_t i, iter = 0; const size_t n = 2; struct rparams p = {1.0, 10.0}; gsl_multiroot_function f = {&rosenbrock_f, n, &p}; double x_init[2] = {-10.0, -5.0}; gsl_vector *x = gsl_vector_alloc (n); gsl_vector_set (x, 0, x_init[0]); gsl_vector_set (x, 1, x_init[1]); T = gsl_multiroot_fsolver_hybrids; s = gsl_multiroot_fsolver_alloc (T, 2); gsl_multiroot_fsolver_set (s, &f, x); print_state (iter, s); do { iter++; status = gsl_multiroot_fsolver_iterate (s); print_state (iter, s); if (status) /* check if solver is stuck */ break; status = gsl_multiroot_test_residual (s->f, 1e-7); } while (status == GSL_CONTINUE && iter < 1000); printf ("status = %s\n", gsl_strerror (status)); gsl_multiroot_fsolver_free (s); gsl_vector_free (x); return 0; } Note that it is important to check the return status of each solver step, in case the algorithm becomes stuck. If an error condition is detected, indicating that the algorithm cannot proceed, then the error can be reported to the user, a new starting point chosen or a different algorithm used. The intermediate state of the solution is displayed by the following function. The solver state contains the vector `s->x' which is the current position, and the vector `s->f' with corresponding function values. int print_state (size_t iter, gsl_multiroot_fsolver * s) { printf ("iter = %3u x = % .3f % .3f " "f(x) = % .3e % .3e\n", iter, gsl_vector_get (s->x, 0), gsl_vector_get (s->x, 1), gsl_vector_get (s->f, 0), gsl_vector_get (s->f, 1)); } Here are the results of running the program. The algorithm is started at (-10,-5) far from the solution. Since the solution is hidden in a narrow valley the earliest steps follow the gradient of the function downhill, in an attempt to reduce the large value of the residual. Once the root has been approximately located, on iteration 8, the Newton behavior takes over and convergence is very rapid. iter = 0 x = -10.000 -5.000 f(x) = 1.100e+01 -1.050e+03 iter = 1 x = -10.000 -5.000 f(x) = 1.100e+01 -1.050e+03 iter = 2 x = -3.976 24.827 f(x) = 4.976e+00 9.020e+01 iter = 3 x = -3.976 24.827 f(x) = 4.976e+00 9.020e+01 iter = 4 x = -3.976 24.827 f(x) = 4.976e+00 9.020e+01 iter = 5 x = -1.274 -5.680 f(x) = 2.274e+00 -7.302e+01 iter = 6 x = -1.274 -5.680 f(x) = 2.274e+00 -7.302e+01 iter = 7 x = 0.249 0.298 f(x) = 7.511e-01 2.359e+00 iter = 8 x = 0.249 0.298 f(x) = 7.511e-01 2.359e+00 iter = 9 x = 1.000 0.878 f(x) = 1.268e-10 -1.218e+00 iter = 10 x = 1.000 0.989 f(x) = 1.124e-11 -1.080e-01 iter = 11 x = 1.000 1.000 f(x) = 0.000e+00 0.000e+00 status = success Note that the algorithm does not update the location on every iteration. Some iterations are used to adjust the trust-region parameter, after trying a step which was found to be divergent, or to recompute the Jacobian, when poor convergence behavior is detected. The next example program adds derivative information, in order to accelerate the solution. There are two derivative functions `rosenbrock_df' and `rosenbrock_fdf'. The latter computes both the function and its derivative simultaneously. This allows the optimization of any common terms. For simplicity we substitute calls to the separate `f' and `df' functions at this point in the code below. int rosenbrock_df (const gsl_vector * x, void *params, gsl_matrix * J) { const double a = ((struct rparams *) params)->a; const double b = ((struct rparams *) params)->b; const double x0 = gsl_vector_get (x, 0); const double df00 = -a; const double df01 = 0; const double df10 = -2 * b * x0; const double df11 = b; gsl_matrix_set (J, 0, 0, df00); gsl_matrix_set (J, 0, 1, df01); gsl_matrix_set (J, 1, 0, df10); gsl_matrix_set (J, 1, 1, df11); return GSL_SUCCESS; } int rosenbrock_fdf (const gsl_vector * x, void *params, gsl_vector * f, gsl_matrix * J) { rosenbrock_f (x, params, f); rosenbrock_df (x, params, J); return GSL_SUCCESS; } The main program now makes calls to the corresponding `fdfsolver' versions of the functions, int main (void) { const gsl_multiroot_fdfsolver_type *T; gsl_multiroot_fdfsolver *s; int status; size_t i, iter = 0; const size_t n = 2; struct rparams p = {1.0, 10.0}; gsl_multiroot_function_fdf f = {&rosenbrock_f, &rosenbrock_df, &rosenbrock_fdf, n, &p}; double x_init[2] = {-10.0, -5.0}; gsl_vector *x = gsl_vector_alloc (n); gsl_vector_set (x, 0, x_init[0]); gsl_vector_set (x, 1, x_init[1]); T = gsl_multiroot_fdfsolver_gnewton; s = gsl_multiroot_fdfsolver_alloc (T, n); gsl_multiroot_fdfsolver_set (s, &f, x); print_state (iter, s); do { iter++; status = gsl_multiroot_fdfsolver_iterate (s); print_state (iter, s); if (status) break; status = gsl_multiroot_test_residual (s->f, 1e-7); } while (status == GSL_CONTINUE && iter < 1000); printf ("status = %s\n", gsl_strerror (status)); gsl_multiroot_fdfsolver_free (s); gsl_vector_free (x); return 0; } The addition of derivative information to the `hybrids' solver does not make any significant difference to its behavior, since it able to approximate the Jacobian numerically with sufficient accuracy. To illustrate the behavior of a different derivative solver we switch to `gnewton'. This is a traditional Newton solver with the constraint that it scales back its step if the full step would lead "uphill". Here is the output for the `gnewton' algorithm, iter = 0 x = -10.000 -5.000 f(x) = 1.100e+01 -1.050e+03 iter = 1 x = -4.231 -65.317 f(x) = 5.231e+00 -8.321e+02 iter = 2 x = 1.000 -26.358 f(x) = -8.882e-16 -2.736e+02 iter = 3 x = 1.000 1.000 f(x) = -2.220e-16 -4.441e-15 status = success The convergence is much more rapid, but takes a wide excursion out to the point (-4.23,-65.3). This could cause the algorithm to go astray in a realistic application. The hybrid algorithm follows the downhill path to the solution more reliably.  File: gsl-ref.info, Node: References and Further Reading for Multidimensional Root Finding, Prev: Example programs for Multidimensional Root finding, Up: Multidimensional Root-Finding 34.9 References and Further Reading =================================== The original version of the Hybrid method is described in the following articles by Powell, M.J.D. Powell, "A Hybrid Method for Nonlinear Equations" (Chap 6, p 87-114) and "A Fortran Subroutine for Solving systems of Nonlinear Algebraic Equations" (Chap 7, p 115-161), in `Numerical Methods for Nonlinear Algebraic Equations', P. Rabinowitz, editor. Gordon and Breach, 1970. The following papers are also relevant to the algorithms described in this section, J.J. More', M.Y. Cosnard, "Numerical Solution of Nonlinear Equations", `ACM Transactions on Mathematical Software', Vol 5, No 1, (1979), p 64-85 C.G. Broyden, "A Class of Methods for Solving Nonlinear Simultaneous Equations", `Mathematics of Computation', Vol 19 (1965), p 577-593 J.J. More', B.S. Garbow, K.E. Hillstrom, "Testing Unconstrained Optimization Software", ACM Transactions on Mathematical Software, Vol 7, No 1 (1981), p 17-41  File: gsl-ref.info, Node: Multidimensional Minimization, Next: Least-Squares Fitting, Prev: Multidimensional Root-Finding, Up: Top 35 Multidimensional Minimization ******************************** This chapter describes routines for finding minima of arbitrary multidimensional 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, while providing 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 minimization algorithms can be used to maximize a function by inverting its sign. The header file `gsl_multimin.h' contains prototypes for the minimization functions and related declarations. * Menu: * Multimin Overview:: * Multimin Caveats:: * Initializing the Multidimensional Minimizer:: * Providing a function to minimize:: * Multimin Iteration:: * Multimin Stopping Criteria:: * Multimin Algorithms:: * Multimin Examples:: * Multimin References and Further Reading::  File: gsl-ref.info, Node: Multimin Overview, Next: Multimin Caveats, Up: Multidimensional Minimization 35.1 Overview ============= The problem of multidimensional minimization requires finding a point x such that the scalar function, f(x_1, ..., x_n) takes a value which is lower than at any neighboring point. For smooth functions the gradient g = \nabla f vanishes at the minimum. In general there are no bracketing methods available for the minimization of n-dimensional functions. The algorithms proceed from an initial guess using a search algorithm which attempts to move in a downhill direction. Algorithms making use of the gradient of the function perform a one-dimensional line minimisation along this direction until the lowest point is found to a suitable tolerance. The search direction is then updated with local information from the function and its derivatives, and the whole process repeated until the true n-dimensional minimum is found. The Nelder-Mead Simplex algorithm applies a different strategy. It maintains n+1 trial parameter vectors as the vertices of a n-dimensional simplex. In each iteration step it tries to improve the worst vertex by a simple geometrical transformation until the size of the simplex falls below a given tolerance. Both types of algorithms use a standard framework. 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 minimizer state, S, for algorithm T * update S using the iteration T * test S for convergence, and repeat iteration if necessary Each iteration step consists either of an improvement to the line-minimisation in the current direction or an update to the search direction itself. The state for the minimizers is held in a `gsl_multimin_fdfminimizer' struct or a `gsl_multimin_fminimizer' struct.  File: gsl-ref.info, Node: Multimin Caveats, Next: Initializing the Multidimensional Minimizer, Prev: Multimin Overview, Up: Multidimensional Minimization 35.2 Caveats ============ Note that the minimization algorithms can only search for one local minimum at a time. When there are several local 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 local minimum in an area where there is more than one. It is also important to note that the minimization algorithms find local minima; there is no way to determine whether a minimum is a global minimum of the function in question.  File: gsl-ref.info, Node: Initializing the Multidimensional Minimizer, Next: Providing a function to minimize, Prev: Multimin Caveats, Up: Multidimensional Minimization 35.3 Initializing the Multidimensional Minimizer ================================================ The following function initializes a multidimensional minimizer. The minimizer itself depends only on the dimension of the problem and the algorithm and can be reused for different problems. -- Function: gsl_multimin_fdfminimizer * gsl_multimin_fdfminimizer_alloc (const gsl_multimin_fdfminimizer_type * T, size_t N) -- Function: gsl_multimin_fminimizer * gsl_multimin_fminimizer_alloc (const gsl_multimin_fminimizer_type * T, size_t N) This function returns a pointer to a newly allocated instance of a minimizer of type T for an N-dimension function. 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_multimin_fdfminimizer_set (gsl_multimin_fdfminimizer * S, gsl_multimin_function_fdf * FDF, const gsl_vector * X, double STEP_SIZE, double TOL) This function initializes the minimizer S to minimize the function FDF starting from the initial point X. The size of the first trial step is given by STEP_SIZE. The accuracy of the line minimization is specified by TOL. The precise meaning of this parameter depends on the method used. Typically the line minimization is considered successful if the gradient of the function g is orthogonal to the current search direction p to a relative accuracy of TOL, where dot(p,g) < tol |p| |g|. A TOL value of 0.1 is suitable for most purposes, since line minimization only needs to be carried out approximately. Note that setting TOL to zero will force the use of "exact" line-searches, which are extremely expensive. -- Function: int gsl_multimin_fminimizer_set (gsl_multimin_fminimizer * S, gsl_multimin_function * F, const gsl_vector * X, const gsl_vector * STEP_SIZE) This function initializes the minimizer S to minimize the function F, starting from the initial point X. The size of the initial trial steps is given in vector STEP_SIZE. The precise meaning of this parameter depends on the method used. -- Function: void gsl_multimin_fdfminimizer_free (gsl_multimin_fdfminimizer * S) -- Function: void gsl_multimin_fminimizer_free (gsl_multimin_fminimizer * S) This function frees all the memory associated with the minimizer S. -- Function: const char * gsl_multimin_fdfminimizer_name (const gsl_multimin_fdfminimizer * S) -- Function: const char * gsl_multimin_fminimizer_name (const gsl_multimin_fminimizer * S) This function returns a pointer to the name of the minimizer. For example, printf ("s is a '%s' minimizer\n", gsl_multimin_fdfminimizer_name (s)); would print something like `s is a 'conjugate_pr' minimizer'.  File: gsl-ref.info, Node: Providing a function to minimize, Next: Multimin Iteration, Prev: Initializing the Multidimensional Minimizer, Up: Multidimensional Minimization 35.4 Providing a function to minimize ===================================== You must provide a parametric function of n variables for the minimizers to operate on. You may also need to provide a routine which calculates the gradient of the function and a third routine which calculates both the function value and the gradient together. In order to allow for general parameters the functions are defined by the following data types: -- Data Type: gsl_multimin_function_fdf This data type defines a general function of n variables with parameters and the corresponding gradient vector of derivatives, `double (* f) (const gsl_vector * X, void * PARAMS)' this function should return the result f(x,params) for argument X and parameters PARAMS. `void (* df) (const gsl_vector * X, void * PARAMS, gsl_vector * G)' this function should store the N-dimensional gradient g_i = d f(x,params) / d x_i in the vector G for argument X and parameters PARAMS, returning an appropriate error code if the function cannot be computed. `void (* fdf) (const gsl_vector * X, void * PARAMS, double * f, gsl_vector * G)' This function should set the values of the F and G as above, for arguments X and parameters PARAMS. This function provides an optimization of the separate functions for f(x) and g(x)--it is always faster to compute the function and its derivative at the same time. `size_t n' the dimension of the system, i.e. the number of components of the vectors X. `void * params' a pointer to the parameters of the function. -- Data Type: gsl_multimin_function This data type defines a general function of n variables with parameters, `double (* f) (const gsl_vector * X, void * PARAMS)' this function should return the result f(x,params) for argument X and parameters PARAMS. `size_t n' the dimension of the system, i.e. the number of components of the vectors X. `void * params' a pointer to the parameters of the function. The following example function defines a simple paraboloid with two parameters, /* Paraboloid centered on (dp[0],dp[1]) */ double my_f (const gsl_vector *v, void *params) { double x, y; double *dp = (double *)params; x = gsl_vector_get(v, 0); y = gsl_vector_get(v, 1); return 10.0 * (x - dp[0]) * (x - dp[0]) + 20.0 * (y - dp[1]) * (y - dp[1]) + 30.0; } /* The gradient of f, df = (df/dx, df/dy). */ void my_df (const gsl_vector *v, void *params, gsl_vector *df) { double x, y; double *dp = (double *)params; x = gsl_vector_get(v, 0); y = gsl_vector_get(v, 1); gsl_vector_set(df, 0, 20.0 * (x - dp[0])); gsl_vector_set(df, 1, 40.0 * (y - dp[1])); } /* Compute both f and df together. */ void my_fdf (const gsl_vector *x, void *params, double *f, gsl_vector *df) { *f = my_f(x, params); my_df(x, params, df); } The function can be initialized using the following code, gsl_multimin_function_fdf my_func; double p[2] = { 1.0, 2.0 }; /* center at (1,2) */ my_func.f = &my_f; my_func.df = &my_df; my_func.fdf = &my_fdf; my_func.n = 2; my_func.params = (void *)p;  File: gsl-ref.info, Node: Multimin Iteration, Next: Multimin Stopping Criteria, Prev: Providing a function to minimize, Up: Multidimensional Minimization 35.5 Iteration ============== The following function drives the iteration of each algorithm. The function performs one iteration to update the state of the minimizer. The same function works for all minimizers so that different methods can be substituted at runtime without modifications to the code. -- Function: int gsl_multimin_fdfminimizer_iterate (gsl_multimin_fdfminimizer * S) -- Function: int gsl_multimin_fminimizer_iterate (gsl_multimin_fminimizer * S) These functions perform a single iteration of the minimizer S. If the iteration encounters an unexpected problem then an error code will be returned. The minimizer maintains a current best estimate of the minimum at all times. This information can be accessed with the following auxiliary functions, -- Function: gsl_vector * gsl_multimin_fdfminimizer_x (const gsl_multimin_fdfminimizer * S) -- Function: gsl_vector * gsl_multimin_fminimizer_x (const gsl_multimin_fminimizer * S) -- Function: double gsl_multimin_fdfminimizer_minimum (const gsl_multimin_fdfminimizer * S) -- Function: double gsl_multimin_fminimizer_minimum (const gsl_multimin_fminimizer * S) -- Function: gsl_vector * gsl_multimin_fdfminimizer_gradient (const gsl_multimin_fdfminimizer * S) -- Function: double gsl_multimin_fminimizer_size (const gsl_multimin_fminimizer * S) These functions return the current best estimate of the location of the minimum, the value of the function at that point, its gradient, and minimizer specific characteristic size for the minimizer S. -- Function: int gsl_multimin_fdfminimizer_restart (gsl_multimin_fdfminimizer * S) This function resets the minimizer S to use the current point as a new starting point.  File: gsl-ref.info, Node: Multimin Stopping Criteria, Next: Multimin Algorithms, Prev: Multimin Iteration, Up: Multidimensional Minimization 35.6 Stopping Criteria ====================== 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 functions below allow the user to test the precision of the current result. -- Function: int gsl_multimin_test_gradient (const gsl_vector * G, double EPSABS) This function tests the norm of the gradient G against the absolute tolerance EPSABS. The gradient of a multidimensional function goes to zero at a minimum. The test returns `GSL_SUCCESS' if the following condition is achieved, |g| < epsabs and returns `GSL_CONTINUE' otherwise. A suitable choice of EPSABS can be made from the desired accuracy in the function for small variations in x. The relationship between these quantities is given by \delta f = g \delta x. -- Function: int gsl_multimin_test_size (const double SIZE, double EPSABS) This function tests the minimizer specific characteristic size (if applicable to the used minimizer) against absolute tolerance EPSABS. The test returns `GSL_SUCCESS' if the size is smaller than tolerance, otherwise `GSL_CONTINUE' is returned.  File: gsl-ref.info, Node: Multimin Algorithms, Next: Multimin Examples, Prev: Multimin Stopping Criteria, Up: Multidimensional Minimization 35.7 Algorithms =============== There are several minimization methods available. The best choice of algorithm depends on the problem. All of the algorithms use the value of the function and its gradient at each evaluation point, except for the Simplex algorithm which uses function values only. -- Minimizer: gsl_multimin_fdfminimizer_conjugate_fr This is the Fletcher-Reeves conjugate gradient algorithm. The conjugate gradient algorithm proceeds as a succession of line minimizations. The sequence of search directions is used to build up an approximation to the curvature of the function in the neighborhood of the minimum. An initial search direction P is chosen using the gradient, and line minimization is carried out in that direction. The accuracy of the line minimization is specified by the parameter TOL. The minimum along this line occurs when the function gradient G and the search direction P are orthogonal. The line minimization terminates when dot(p,g) < tol |p| |g|. The search direction is updated using the Fletcher-Reeves formula p' = g' - \beta g where \beta=-|g'|^2/|g|^2, and the line minimization is then repeated for the new search direction. -- Minimizer: gsl_multimin_fdfminimizer_conjugate_pr This is the Polak-Ribiere conjugate gradient algorithm. It is similar to the Fletcher-Reeves method, differing only in the choice of the coefficient \beta. Both methods work well when the evaluation point is close enough to the minimum of the objective function that it is well approximated by a quadratic hypersurface. -- Minimizer: gsl_multimin_fdfminimizer_vector_bfgs2 -- Minimizer: gsl_multimin_fdfminimizer_vector_bfgs These methods use the vector Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm. This is a quasi-Newton method which builds up an approximation to the second derivatives of the function f using the difference between successive gradient vectors. By combining the first and second derivatives the algorithm is able to take Newton-type steps towards the function minimum, assuming quadratic behavior in that region. The `bfgs2' version of this minimizer is the most efficient version available, and is a faithful implementation of the line minimization scheme described in Fletcher's `Practical Methods of Optimization', Algorithms 2.6.2 and 2.6.4. It supercedes the original `bfgs' routine and requires substantially fewer function and gradient evaluations. The user-supplied tolerance TOL corresponds to the parameter \sigma used by Fletcher. A value of 0.1 is recommended for typical use (larger values correspond to less accurate line searches). -- Minimizer: gsl_multimin_fdfminimizer_steepest_descent The steepest descent algorithm follows the downhill gradient of the function at each step. When a downhill step is successful the step-size is increased by a factor of two. If the downhill step leads to a higher function value then the algorithm backtracks and the step size is decreased using the parameter TOL. A suitable value of TOL for most applications is 0.1. The steepest descent method is inefficient and is included only for demonstration purposes. -- Minimizer: gsl_multimin_fminimizer_nmsimplex This is the Simplex algorithm of Nelder and Mead. It constructs n vectors p_i from the starting vector X and the vector STEP_SIZE as follows: p_0 = (x_0, x_1, ... , x_n) p_1 = (x_0 + step_size_0, x_1, ... , x_n) p_2 = (x_0, x_1 + step_size_1, ... , x_n) ... = ... p_n = (x_0, x_1, ... , x_n+step_size_n) These vectors form the n+1 vertices of a simplex in n dimensions. On each iteration the algorithm tries to improve the parameter vector p_i corresponding to the highest function value by simple geometrical transformations. These are reflection, reflection followed by expansion, contraction and multiple contraction. Using these transformations the simplex moves through the parameter space towards the minimum, where it contracts itself. After each iteration, the best vertex is returned. Note, that due to the nature of the algorithm not every step improves the current best parameter vector. Usually several iterations are required. The routine calculates the minimizer specific characteristic size as the average distance from the geometrical center of the simplex to all its vertices. This size can be used as a stopping criteria, as the simplex contracts itself near the minimum. The size is returned by the function `gsl_multimin_fminimizer_size'.  File: gsl-ref.info, Node: Multimin Examples, Next: Multimin References and Further Reading, Prev: Multimin Algorithms, Up: Multidimensional Minimization 35.8 Examples ============= This example program finds the minimum of the paraboloid function defined earlier. The location of the minimum is offset from the origin in x and y, and the function value at the minimum is non-zero. The main program is given below, it requires the example function given earlier in this chapter. int main (void) { size_t iter = 0; int status; const gsl_multimin_fdfminimizer_type *T; gsl_multimin_fdfminimizer *s; /* Position of the minimum (1,2). */ double par[2] = { 1.0, 2.0 }; gsl_vector *x; gsl_multimin_function_fdf my_func; my_func.f = &my_f; my_func.df = &my_df; my_func.fdf = &my_fdf; my_func.n = 2; my_func.params = ∥ /* Starting point, x = (5,7) */ x = gsl_vector_alloc (2); gsl_vector_set (x, 0, 5.0); gsl_vector_set (x, 1, 7.0); T = gsl_multimin_fdfminimizer_conjugate_fr; s = gsl_multimin_fdfminimizer_alloc (T, 2); gsl_multimin_fdfminimizer_set (s, &my_func, x, 0.01, 1e-4); do { iter++; status = gsl_multimin_fdfminimizer_iterate (s); if (status) break; status = gsl_multimin_test_gradient (s->gradient, 1e-3); if (status == GSL_SUCCESS) printf ("Minimum found at:\n"); printf ("%5d %.5f %.5f %10.5f\n", iter, gsl_vector_get (s->x, 0), gsl_vector_get (s->x, 1), s->f); } while (status == GSL_CONTINUE && iter < 100); gsl_multimin_fdfminimizer_free (s); gsl_vector_free (x); return 0; } The initial step-size is chosen as 0.01, a conservative estimate in this case, and the line minimization parameter is set at 0.0001. The program terminates when the norm of the gradient has been reduced below 0.001. The output of the program is shown below, x y f 1 4.99629 6.99072 687.84780 2 4.98886 6.97215 683.55456 3 4.97400 6.93501 675.01278 4 4.94429 6.86073 658.10798 5 4.88487 6.71217 625.01340 6 4.76602 6.41506 561.68440 7 4.52833 5.82083 446.46694 8 4.05295 4.63238 261.79422 9 3.10219 2.25548 75.49762 10 2.85185 1.62963 67.03704 11 2.19088 1.76182 45.31640 12 0.86892 2.02622 30.18555 Minimum found at: 13 1.00000 2.00000 30.00000 Note that the algorithm gradually increases the step size as it successfully moves downhill, as can be seen by plotting the successive points. The conjugate gradient algorithm finds the minimum on its second direction because the function is purely quadratic. Additional iterations would be needed for a more complicated function. Here is another example using the Nelder-Mead Simplex algorithm to minimize the same example object function, as above. int main(void) { size_t np = 2; double par[2] = {1.0, 2.0}; const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex; gsl_multimin_fminimizer *s = NULL; gsl_vector *ss, *x; gsl_multimin_function minex_func; size_t iter = 0, i; int status; double size; /* Initial vertex size vector */ ss = gsl_vector_alloc (np); /* Set all step sizes to 1 */ gsl_vector_set_all (ss, 1.0); /* Starting point */ x = gsl_vector_alloc (np); gsl_vector_set (x, 0, 5.0); gsl_vector_set (x, 1, 7.0); /* Initialize method and iterate */ minex_func.f = &my_f; minex_func.n = np; minex_func.params = (void *)∥ s = gsl_multimin_fminimizer_alloc (T, np); gsl_multimin_fminimizer_set (s, &minex_func, x, ss); do { iter++; status = gsl_multimin_fminimizer_iterate(s); if (status) break; size = gsl_multimin_fminimizer_size (s); status = gsl_multimin_test_size (size, 1e-2); if (status == GSL_SUCCESS) { printf ("converged to minimum at\n"); } printf ("%5d ", iter); for (i = 0; i < np; i++) { printf ("%10.3e ", gsl_vector_get (s->x, i)); } printf ("f() = %7.3f size = %.3f\n", s->fval, size); } while (status == GSL_CONTINUE && iter < 100); gsl_vector_free(x); gsl_vector_free(ss); gsl_multimin_fminimizer_free (s); return status; } The minimum search stops when the Simplex size drops to 0.01. The output is shown below. 1 6.500e+00 5.000e+00 f() = 512.500 size = 1.082 2 5.250e+00 4.000e+00 f() = 290.625 size = 1.372 3 5.250e+00 4.000e+00 f() = 290.625 size = 1.372 4 5.500e+00 1.000e+00 f() = 252.500 size = 1.372 5 2.625e+00 3.500e+00 f() = 101.406 size = 1.823 6 3.469e+00 1.375e+00 f() = 98.760 size = 1.526 7 1.820e+00 3.156e+00 f() = 63.467 size = 1.105 8 1.820e+00 3.156e+00 f() = 63.467 size = 1.105 9 1.016e+00 2.812e+00 f() = 43.206 size = 1.105 10 2.041e+00 2.008e+00 f() = 40.838 size = 0.645 11 1.236e+00 1.664e+00 f() = 32.816 size = 0.645 12 1.236e+00 1.664e+00 f() = 32.816 size = 0.447 13 5.225e-01 1.980e+00 f() = 32.288 size = 0.447 14 1.103e+00 2.073e+00 f() = 30.214 size = 0.345 15 1.103e+00 2.073e+00 f() = 30.214 size = 0.264 16 1.103e+00 2.073e+00 f() = 30.214 size = 0.160 17 9.864e-01 1.934e+00 f() = 30.090 size = 0.132 18 9.190e-01 1.987e+00 f() = 30.069 size = 0.092 19 1.028e+00 2.017e+00 f() = 30.013 size = 0.056 20 1.028e+00 2.017e+00 f() = 30.013 size = 0.046 21 1.028e+00 2.017e+00 f() = 30.013 size = 0.033 22 9.874e-01 1.985e+00 f() = 30.006 size = 0.028 23 9.846e-01 1.995e+00 f() = 30.003 size = 0.023 24 1.007e+00 2.003e+00 f() = 30.001 size = 0.012 converged to minimum at 25 1.007e+00 2.003e+00 f() = 30.001 size = 0.010 The simplex size first increases, while the simplex moves towards the minimum. After a while the size begins to decrease as the simplex contracts around the minimum.  File: gsl-ref.info, Node: Multimin References and Further Reading, Prev: Multimin Examples, Up: Multidimensional Minimization 35.9 References and Further Reading =================================== The conjugate gradient and BFGS methods are described in detail in the following book, R. Fletcher, `Practical Methods of Optimization (Second Edition)' Wiley (1987), ISBN 0471915475. A brief description of multidimensional minimization algorithms and more recent references can be found in, C.W. Ueberhuber, `Numerical Computation (Volume 2)', Chapter 14, Section 4.4 "Minimization Methods", p. 325-335, Springer (1997), ISBN 3-540-62057-5. The simplex algorithm is described in the following paper, J.A. Nelder and R. Mead, `A simplex method for function minimization', Computer Journal vol. 7 (1965), 308-315.  File: gsl-ref.info, Node: Least-Squares Fitting, Next: Nonlinear Least-Squares Fitting, Prev: Multidimensional Minimization, Up: Top 36 Least-Squares Fitting ************************ This chapter describes routines for performing least squares fits to experimental data using linear combinations of functions. The data may be weighted or unweighted, i.e. with known or unknown errors. For weighted data the functions compute the best fit parameters and their associated covariance matrix. For unweighted data the covariance matrix is estimated from the scatter of the points, giving a variance-covariance matrix. The functions are divided into separate versions for simple one- or two-parameter regression and multiple-parameter fits. The functions are declared in the header file `gsl_fit.h'. * Menu: * Fitting Overview:: * Linear regression:: * Linear fitting without a constant term:: * Multi-parameter fitting:: * Fitting Examples:: * Fitting References and Further Reading::  File: gsl-ref.info, Node: Fitting Overview, Next: Linear regression, Up: Least-Squares Fitting 36.1 Overview ============= Least-squares fits are found by minimizing \chi^2 (chi-squared), the weighted sum of squared residuals over n experimental datapoints (x_i, y_i) for the model Y(c,x), \chi^2 = \sum_i w_i (y_i - Y(c, x_i))^2 The p parameters of the model are c = {c_0, c_1, ...}. The weight factors w_i are given by w_i = 1/\sigma_i^2, where \sigma_i is the experimental error on the data-point y_i. The errors are assumed to be gaussian and uncorrelated. For unweighted data the chi-squared sum is computed without any weight factors. The fitting routines return the best-fit parameters c and their p \times p covariance matrix. The covariance matrix measures the statistical errors on the best-fit parameters resulting from the errors on the data, \sigma_i, and is defined as C_{ab} = <\delta c_a \delta c_b> where < > denotes an average over the gaussian error distributions of the underlying datapoints. The covariance matrix is calculated by error propagation from the data errors \sigma_i. The change in a fitted parameter \delta c_a caused by a small change in the data \delta y_i is given by \delta c_a = \sum_i (dc_a/dy_i) \delta y_i allowing the covariance matrix to be written in terms of the errors on the data, C_{ab} = \sum_{i,j} (dc_a/dy_i) (dc_b/dy_j) <\delta y_i \delta y_j> For uncorrelated data the fluctuations of the underlying datapoints satisfy <\delta y_i \delta y_j> = \sigma_i^2 \delta_{ij}, giving a corresponding parameter covariance matrix of C_{ab} = \sum_i (1/w_i) (dc_a/dy_i) (dc_b/dy_i) When computing the covariance matrix for unweighted data, i.e. data with unknown errors, the weight factors w_i in this sum are replaced by the single estimate w = 1/\sigma^2, where \sigma^2 is the computed variance of the residuals about the best-fit model, \sigma^2 = \sum (y_i - Y(c,x_i))^2 / (n-p). This is referred to as the "variance-covariance matrix". The standard deviations of the best-fit parameters are given by the square root of the corresponding diagonal elements of the covariance matrix, \sigma_{c_a} = \sqrt{C_{aa}}.  File: gsl-ref.info, Node: Linear regression, Next: Linear fitting without a constant term, Prev: Fitting Overview, Up: Least-Squares Fitting 36.2 Linear regression ====================== The functions described in this section can be used to perform least-squares fits to a straight line model, Y(c,x) = c_0 + c_1 x. -- Function: int gsl_fit_linear (const double * X, const size_t XSTRIDE, const double * Y, const size_t YSTRIDE, size_t N, double * C0, double * C1, double * COV00, double * COV01, double * COV11, double * SUMSQ) This function computes the best-fit linear regression coefficients (C0,C1) of the model Y = c_0 + c_1 X for the dataset (X, Y), two vectors of length N with strides XSTRIDE and YSTRIDE. The errors on Y are assumed unknown so the variance-covariance matrix for the parameters (C0, C1) is estimated from the scatter of the points around the best-fit line and returned via the parameters (COV00, COV01, COV11). The sum of squares of the residuals from the best-fit line is returned in SUMSQ. -- Function: int gsl_fit_wlinear (const double * X, const size_t XSTRIDE, const double * W, const size_t WSTRIDE, const double * Y, const size_t YSTRIDE, size_t N, double * C0, double * C1, double * COV00, double * COV01, double * COV11, double * CHISQ) This function computes the best-fit linear regression coefficients (C0,C1) of the model Y = c_0 + c_1 X for the weighted dataset (X, Y), two vectors of length N with strides XSTRIDE and YSTRIDE. The vector W, of length N and stride WSTRIDE, specifies the weight of each datapoint. The weight is the reciprocal of the variance for each datapoint in Y. The covariance matrix for the parameters (C0, C1) is computed using the weights and returned via the parameters (COV00, COV01, COV11). The weighted sum of squares of the residuals from the best-fit line, \chi^2, is returned in CHISQ. -- Function: int gsl_fit_linear_est (double X, double C0, double C1, double C00, double C01, double C11, double * Y, double * Y_ERR) This function uses the best-fit linear regression coefficients C0,C1 and their covariance COV00,COV01,COV11 to compute the fitted function Y and its standard deviation Y_ERR for the model Y = c_0 + c_1 X at the point X.