function/************************************************************************* K-NN query: K nearest neighbors INPUT PARAMETERS KDT - KD-tree X - point, array[0..NX-1]. K - number of neighbors to return, K>=1 SelfMatch - whether self-matches are allowed: * if True, nearest neighbor may be the point itself (if it exists in original dataset) * if False, then only points with non-zero distance are returned * if not given, considered True RESULT number of actual neighbors found (either K or N, if K>N). This subroutine performs query and stores its result in the internal structures of the KD-tree. You can use following subroutines to obtain these results: * KDTreeQueryResultsX() to get X-values * KDTreeQueryResultsXY() to get X- and Y-values * KDTreeQueryResultsTags() to get tag values * KDTreeQueryResultsDistances() to get distances -- ALGLIB -- Copyright 28.02.2010 by Bochkanov Sergey *************************************************************************/public static int alglib.kdtreequeryknn(kdtree kdt, double[] x, int k) public static int alglib.kdtreequeryknn( kdtree kdt, double[] x, int k, bool selfmatch)
Examples: [1]
kdtreequeryresultsdistances
function/************************************************************************* Distances from last query INPUT PARAMETERS KDT - KD-tree R - possibly pre-allocated buffer. If X is too small to store result, it is resized. If size(X) is enough to store result, it is left unchanged. OUTPUT PARAMETERS R - filled with distances (in corresponding norm) NOTES 1. points are ordered by distance from the query point (first = closest) 2. if XY is larger than required to store result, only leading part will be overwritten; trailing part will be left unchanged. So if on input XY = [[A,B],[C,D]], and result is [1,2], then on exit we will get XY = [[1,2],[C,D]]. This is done purposely to increase performance; if you want function to resize array according to result size, use function with same name and suffix 'I'. SEE ALSO * KDTreeQueryResultsX() X-values * KDTreeQueryResultsXY() X- and Y-values * KDTreeQueryResultsTags() tag values -- ALGLIB -- Copyright 28.02.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.kdtreequeryresultsdistances( kdtree kdt, ref double[] r)
Examples: [1]
kdtreequeryresultsdistancesi
function/************************************************************************* Distances from last query; 'interactive' variant for languages like Python which support constructs like "R = KDTreeQueryResultsDistancesI(KDT)" and interactive mode of interpreter. This function allocates new array on each call, so it is significantly slower than its 'non-interactive' counterpart, but it is more convenient when you call it from command line. -- ALGLIB -- Copyright 28.02.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.kdtreequeryresultsdistancesi( kdtree kdt, out double[] r)
kdtreequeryresultstags
function/************************************************************************* Tags from last query INPUT PARAMETERS KDT - KD-tree Tags - possibly pre-allocated buffer. If X is too small to store result, it is resized. If size(X) is enough to store result, it is left unchanged. OUTPUT PARAMETERS Tags - filled with tags associated with points, or, when no tags were supplied, with zeros NOTES 1. points are ordered by distance from the query point (first = closest) 2. if XY is larger than required to store result, only leading part will be overwritten; trailing part will be left unchanged. So if on input XY = [[A,B],[C,D]], and result is [1,2], then on exit we will get XY = [[1,2],[C,D]]. This is done purposely to increase performance; if you want function to resize array according to result size, use function with same name and suffix 'I'. SEE ALSO * KDTreeQueryResultsX() X-values * KDTreeQueryResultsXY() X- and Y-values * KDTreeQueryResultsDistances() distances -- ALGLIB -- Copyright 28.02.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.kdtreequeryresultstags( kdtree kdt, ref int[] tags)
Examples: [1]
kdtreequeryresultstagsi
function/************************************************************************* Tags from last query; 'interactive' variant for languages like Python which support constructs like "Tags = KDTreeQueryResultsTagsI(KDT)" and interactive mode of interpreter. This function allocates new array on each call, so it is significantly slower than its 'non-interactive' counterpart, but it is more convenient when you call it from command line. -- ALGLIB -- Copyright 28.02.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.kdtreequeryresultstagsi( kdtree kdt, out int[] tags)
kdtreequeryresultsx
function/************************************************************************* X-values from last query INPUT PARAMETERS KDT - KD-tree X - possibly pre-allocated buffer. If X is too small to store result, it is resized. If size(X) is enough to store result, it is left unchanged. OUTPUT PARAMETERS X - rows are filled with X-values NOTES 1. points are ordered by distance from the query point (first = closest) 2. if XY is larger than required to store result, only leading part will be overwritten; trailing part will be left unchanged. So if on input XY = [[A,B],[C,D]], and result is [1,2], then on exit we will get XY = [[1,2],[C,D]]. This is done purposely to increase performance; if you want function to resize array according to result size, use function with same name and suffix 'I'. SEE ALSO * KDTreeQueryResultsXY() X- and Y-values * KDTreeQueryResultsTags() tag values * KDTreeQueryResultsDistances() distances -- ALGLIB -- Copyright 28.02.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.kdtreequeryresultsx(kdtree kdt, ref double[,] x)
Examples: [1]
kdtreequeryresultsxi
function/************************************************************************* X-values from last query; 'interactive' variant for languages like Python which support constructs like "X = KDTreeQueryResultsXI(KDT)" and interactive mode of interpreter. This function allocates new array on each call, so it is significantly slower than its 'non-interactive' counterpart, but it is more convenient when you call it from command line. -- ALGLIB -- Copyright 28.02.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.kdtreequeryresultsxi( kdtree kdt, out double[,] x)
kdtreequeryresultsxy
function/************************************************************************* X- and Y-values from last query INPUT PARAMETERS KDT - KD-tree XY - possibly pre-allocated buffer. If XY is too small to store result, it is resized. If size(XY) is enough to store result, it is left unchanged. OUTPUT PARAMETERS XY - rows are filled with points: first NX columns with X-values, next NY columns - with Y-values. NOTES 1. points are ordered by distance from the query point (first = closest) 2. if XY is larger than required to store result, only leading part will be overwritten; trailing part will be left unchanged. So if on input XY = [[A,B],[C,D]], and result is [1,2], then on exit we will get XY = [[1,2],[C,D]]. This is done purposely to increase performance; if you want function to resize array according to result size, use function with same name and suffix 'I'. SEE ALSO * KDTreeQueryResultsX() X-values * KDTreeQueryResultsTags() tag values * KDTreeQueryResultsDistances() distances -- ALGLIB -- Copyright 28.02.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.kdtreequeryresultsxy( kdtree kdt, ref double[,] xy)
Examples: [1]
kdtreequeryresultsxyi
function/************************************************************************* XY-values from last query; 'interactive' variant for languages like Python which support constructs like "XY = KDTreeQueryResultsXYI(KDT)" and interactive mode of interpreter. This function allocates new array on each call, so it is significantly slower than its 'non-interactive' counterpart, but it is more convenient when you call it from command line. -- ALGLIB -- Copyright 28.02.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.kdtreequeryresultsxyi( kdtree kdt, out double[,] xy)
kdtreequeryrnn
function/************************************************************************* R-NN query: all points within R-sphere centered at X INPUT PARAMETERS KDT - KD-tree X - point, array[0..NX-1]. R - radius of sphere (in corresponding norm), R>0 SelfMatch - whether self-matches are allowed: * if True, nearest neighbor may be the point itself (if it exists in original dataset) * if False, then only points with non-zero distance are returned * if not given, considered True RESULT number of neighbors found, >=0 This subroutine performs query and stores its result in the internal structures of the KD-tree. You can use following subroutines to obtain actual results: * KDTreeQueryResultsX() to get X-values * KDTreeQueryResultsXY() to get X- and Y-values * KDTreeQueryResultsTags() to get tag values * KDTreeQueryResultsDistances() to get distances -- ALGLIB -- Copyright 28.02.2010 by Bochkanov Sergey *************************************************************************/public static int alglib.kdtreequeryrnn(kdtree kdt, double[] x, double r) public static int alglib.kdtreequeryrnn( kdtree kdt, double[] x, double r, bool selfmatch)
Examples: [1]
kdtreeserialize
function/************************************************************************* This function serializes data structure to string. Important properties of s_out: * it contains alphanumeric characters, dots, underscores, minus signs * these symbols are grouped into words, which are separated by spaces and Windows-style (CR+LF) newlines * although serializer uses spaces and CR+LF as separators, you can replace any separator character by arbitrary combination of spaces, tabs, Windows or Unix newlines. It allows flexible reformatting of the string in case you want to include it into text or XML file. But you should not insert separators into the middle of the "words" nor you should change case of letters. * s_out can be freely moved between 32-bit and 64-bit systems, little and big endian machines, and so on. You can serialize structure on 32-bit machine and unserialize it on 64-bit one (or vice versa), or serialize it on SPARC and unserialize on x86. You can also serialize it in C++ version of ALGLIB and unserialize in C# one, and vice versa. *************************************************************************/public static void kdtreeserialize(kdtree obj, out string s_out)
kdtreeunserialize
function/************************************************************************* This function unserializes data structure from string. *************************************************************************/public static void kdtreeunserialize(string s_in, out kdtree obj)
public static int Main(string[] args) { double[,] a = new double[,]{{0,0},{0,1},{1,0},{1,1}}; int nx = 2; int ny = 0; int normtype = 2; alglib.kdtree kdt; double[] x; double[,] r = new double[0,0]; int k; alglib.kdtreebuild(a, nx, ny, normtype, out kdt); x = new double[]{-1,0}; k = alglib.kdtreequeryknn(kdt, x, 1); System.Console.WriteLine("{0}", k); // EXPECTED: 1 alglib.kdtreequeryresultsx(kdt, ref r); System.Console.WriteLine("{0}", alglib.ap.format(r,1)); // EXPECTED: [[0,0]] System.Console.ReadLine(); return 0; }
public static int Main(string[] args) { double[,] a = new double[,]{{0,0},{0,1},{1,0},{1,1}}; int nx = 2; int ny = 0; int normtype = 2; alglib.kdtree kdt0; alglib.kdtree kdt1; string s; double[] x; double[,] r0 = new double[0,0]; double[,] r1 = new double[0,0]; // // Build tree and serialize it // alglib.kdtreebuild(a, nx, ny, normtype, out kdt0); alglib.kdtreeserialize(kdt0, out s); alglib.kdtreeunserialize(s, out kdt1); // // Compare results from KNN queries // x = new double[]{-1,0}; alglib.kdtreequeryknn(kdt0, x, 1); alglib.kdtreequeryresultsx(kdt0, ref r0); alglib.kdtreequeryknn(kdt1, x, 1); alglib.kdtreequeryresultsx(kdt1, ref r1); System.Console.WriteLine("{0}", alglib.ap.format(r0,1)); // EXPECTED: [[0,0]] System.Console.WriteLine("{0}", alglib.ap.format(r1,1)); // EXPECTED: [[0,0]] System.Console.ReadLine(); return 0; }
nleq
subpackagenleqreport
class/************************************************************************* *************************************************************************/public class nleqreport { public int iterationscount; public int nfunc; public int njac; public int terminationtype; }
nleqstate
class/************************************************************************* *************************************************************************/public class nleqstate { }
nleqcreatelm
function/************************************************************************* LEVENBERG-MARQUARDT-LIKE NONLINEAR SOLVER DESCRIPTION: This algorithm solves system of nonlinear equations F[0](x[0], ..., x[n-1]) = 0 F[1](x[0], ..., x[n-1]) = 0 ... F[M-1](x[0], ..., x[n-1]) = 0 with M/N do not necessarily coincide. Algorithm converges quadratically under following conditions: * the solution set XS is nonempty * for some xs in XS there exist such neighbourhood N(xs) that: * vector function F(x) and its Jacobian J(x) are continuously differentiable on N * ||F(x)|| provides local error bound on N, i.e. there exists such c1, that ||F(x)||>c1*distance(x,XS) Note that these conditions are much more weaker than usual non-singularity conditions. For example, algorithm will converge for any affine function F (whether its Jacobian singular or not). REQUIREMENTS: Algorithm will request following information during its operation: * function vector F[] and Jacobian matrix at given point X * value of merit function f(x)=F[0]^2(x)+...+F[M-1]^2(x) at given point X USAGE: 1. User initializes algorithm state with NLEQCreateLM() call 2. User tunes solver parameters with NLEQSetCond(), NLEQSetStpMax() and other functions 3. User calls NLEQSolve() function which takes algorithm state and pointers (delegates, etc.) to callback functions which calculate merit function value and Jacobian. 4. User calls NLEQResults() to get solution 5. Optionally, user may call NLEQRestartFrom() to solve another problem with same parameters (N/M) but another starting point and/or another function vector. NLEQRestartFrom() allows to reuse already initialized structure. INPUT PARAMETERS: N - space dimension, N>1: * if provided, only leading N elements of X are used * if not provided, determined automatically from size of X M - system size X - starting point OUTPUT PARAMETERS: State - structure which stores algorithm state NOTES: 1. you may tune stopping conditions with NLEQSetCond() function 2. if target function contains exp() or other fast growing functions, and optimization algorithm makes too large steps which leads to overflow, use NLEQSetStpMax() function to bound algorithm's steps. 3. this algorithm is a slightly modified implementation of the method described in 'Levenberg-Marquardt method for constrained nonlinear equations with strong local convergence properties' by Christian Kanzow Nobuo Yamashita and Masao Fukushima and further developed in 'On the convergence of a New Levenberg-Marquardt Method' by Jin-yan Fan and Ya-Xiang Yuan. -- ALGLIB -- Copyright 20.08.2009 by Bochkanov Sergey *************************************************************************/public static void alglib.nleqcreatelm( int m, double[] x, out nleqstate state) public static void alglib.nleqcreatelm( int n, int m, double[] x, out nleqstate state)
nleqrestartfrom
function/************************************************************************* This subroutine restarts CG algorithm from new point. All optimization parameters are left unchanged. This function allows to solve multiple optimization problems (which must have same number of dimensions) without object reallocation penalty. INPUT PARAMETERS: State - structure used for reverse communication previously allocated with MinCGCreate call. X - new starting point. BndL - new lower bounds BndU - new upper bounds -- ALGLIB -- Copyright 30.07.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.nleqrestartfrom(nleqstate state, double[] x)
nleqresults
function/************************************************************************* NLEQ solver results INPUT PARAMETERS: State - algorithm state. OUTPUT PARAMETERS: X - array[0..N-1], solution Rep - optimization report: * Rep.TerminationType completetion code: * -4 ERROR: algorithm has converged to the stationary point Xf which is local minimum of f=F[0]^2+...+F[m-1]^2, but is not solution of nonlinear system. * 1 sqrt(f)<=EpsF. * 5 MaxIts steps was taken * 7 stopping conditions are too stringent, further improvement is impossible * Rep.IterationsCount contains iterations count * NFEV countains number of function calculations * ActiveConstraints contains number of active constraints -- ALGLIB -- Copyright 20.08.2009 by Bochkanov Sergey *************************************************************************/public static void alglib.nleqresults( nleqstate state, out double[] x, out nleqreport rep)
nleqresultsbuf
function/************************************************************************* NLEQ solver results Buffered implementation of NLEQResults(), which uses pre-allocated buffer to store X[]. If buffer size is too small, it resizes buffer. It is intended to be used in the inner cycles of performance critical algorithms where array reallocation penalty is too large to be ignored. -- ALGLIB -- Copyright 20.08.2009 by Bochkanov Sergey *************************************************************************/public static void alglib.nleqresultsbuf( nleqstate state, ref double[] x, nleqreport rep)
nleqsetcond
function/************************************************************************* This function sets stopping conditions for the nonlinear solver INPUT PARAMETERS: State - structure which stores algorithm state EpsF - >=0 The subroutine finishes its work if on k+1-th iteration the condition ||F||<=EpsF is satisfied MaxIts - maximum number of iterations. If MaxIts=0, the number of iterations is unlimited. Passing EpsF=0 and MaxIts=0 simultaneously will lead to automatic stopping criterion selection (small EpsF). NOTES: -- ALGLIB -- Copyright 20.08.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.nleqsetcond( nleqstate state, double epsf, int maxits)
nleqsetstpmax
function/************************************************************************* This function sets maximum step length INPUT PARAMETERS: State - structure which stores algorithm state StpMax - maximum step length, >=0. Set StpMax to 0.0, if you don't want to limit step length. Use this subroutine when target function contains exp() or other fast growing functions, and algorithm makes too large steps which lead to overflow. This function allows us to reject steps that are too large (and therefore expose us to the possible overflow) without actually calculating function value at the x+stp*d. -- ALGLIB -- Copyright 20.08.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.nleqsetstpmax(nleqstate state, double stpmax)
nleqsetxrep
function/************************************************************************* This function turns on/off reporting. INPUT PARAMETERS: State - structure which stores algorithm state NeedXRep- whether iteration reports are needed or not If NeedXRep is True, algorithm will call rep() callback function if it is provided to NLEQSolve(). -- ALGLIB -- Copyright 20.08.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.nleqsetxrep(nleqstate state, bool needxrep)
nleqsolve
function/************************************************************************* This family of functions is used to launcn iterations of nonlinear solver These functions accept following parameters: state - algorithm state func - callback which calculates function (or merit function) value func at given point x jac - callback which calculates function vector fi[] and Jacobian jac at given point x rep - optional callback which is called after each iteration can be NULL ptr - optional pointer which is passed to func/grad/hess/jac/rep can be NULL -- ALGLIB -- Copyright 20.03.2009 by Bochkanov Sergey *************************************************************************/public static void nleqsolve(nleqstate state, ndimensional_func func, ndimensional_jac jac, ndimensional_rep rep, object obj)
normaldistr
subpackageerrorfunction
function/************************************************************************* Error function The integral is x - 2 | | 2 erf(x) = -------- | exp( - t ) dt. sqrt(pi) | | - 0 For 0 <= |x| < 1, erf(x) = x * P4(x**2)/Q5(x**2); otherwise erf(x) = 1 - erfc(x). ACCURACY: Relative error: arithmetic domain # trials peak rms IEEE 0,1 30000 3.7e-16 1.0e-16 Cephes Math Library Release 2.8: June, 2000 Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier *************************************************************************/public static double alglib.errorfunction(double x)
errorfunctionc
function/************************************************************************* Complementary error function 1 - erf(x) = inf. - 2 | | 2 erfc(x) = -------- | exp( - t ) dt sqrt(pi) | | - x For small x, erfc(x) = 1 - erf(x); otherwise rational approximations are computed. ACCURACY: Relative error: arithmetic domain # trials peak rms IEEE 0,26.6417 30000 5.7e-14 1.5e-14 Cephes Math Library Release 2.8: June, 2000 Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier *************************************************************************/public static double alglib.errorfunctionc(double x)
inverf
function/************************************************************************* Inverse of the error function Cephes Math Library Release 2.8: June, 2000 Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier *************************************************************************/public static double alglib.inverf(double e)
invnormaldistribution
function/************************************************************************* Inverse of Normal distribution function Returns the argument, x, for which the area under the Gaussian probability density function (integrated from minus infinity to x) is equal to y. For small arguments 0 < y < exp(-2), the program computes z = sqrt( -2.0 * log(y) ); then the approximation is x = z - log(z)/z - (1/z) P(1/z) / Q(1/z). There are two rational functions P/Q, one for 0 < y < exp(-32) and the other for y up to exp(-2). For larger arguments, w = y - 0.5, and x/sqrt(2pi) = w + w**3 R(w**2)/S(w**2)). ACCURACY: Relative error: arithmetic domain # trials peak rms IEEE 0.125, 1 20000 7.2e-16 1.3e-16 IEEE 3e-308, 0.135 50000 4.6e-16 9.8e-17 Cephes Math Library Release 2.8: June, 2000 Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier *************************************************************************/public static double alglib.invnormaldistribution(double y0)
normaldistribution
function/************************************************************************* Normal distribution function Returns the area under the Gaussian probability density function, integrated from minus infinity to x: x - 1 | | 2 ndtr(x) = --------- | exp( - t /2 ) dt sqrt(2pi) | | - -inf. = ( 1 + erf(z) ) / 2 = erfc(z) / 2 where z = x/sqrt(2). Computation is via the functions erf and erfc. ACCURACY: Relative error: arithmetic domain # trials peak rms IEEE -13,0 30000 3.4e-14 6.7e-15 Cephes Math Library Release 2.8: June, 2000 Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier *************************************************************************/public static double alglib.normaldistribution(double x)
normestimator
subpackagenormestimatorstate
class/************************************************************************* This object stores state of the iterative norm estimation algorithm. You should use ALGLIB functions to work with this object. *************************************************************************/public class normestimatorstate { }
normestimatorcreate
function/************************************************************************* This procedure initializes matrix norm estimator. USAGE: 1. User initializes algorithm state with NormEstimatorCreate() call 2. User calls NormEstimatorEstimateSparse() (or NormEstimatorIteration()) 3. User calls NormEstimatorResults() to get solution. INPUT PARAMETERS: M - number of rows in the matrix being estimated, M>0 N - number of columns in the matrix being estimated, N>0 NStart - number of random starting vectors recommended value - at least 5. NIts - number of iterations to do with best starting vector recommended value - at least 5. OUTPUT PARAMETERS: State - structure which stores algorithm state NOTE: this algorithm is effectively deterministic, i.e. it always returns same result when repeatedly called for the same matrix. In fact, algorithm uses randomized starting vectors, but internal random numbers generator always generates same sequence of the random values (it is a feature, not bug). Algorithm can be made non-deterministic with NormEstimatorSetSeed(0) call. -- ALGLIB -- Copyright 06.12.2011 by Bochkanov Sergey *************************************************************************/public static void alglib.normestimatorcreate( int m, int n, int nstart, int nits, out normestimatorstate state)
normestimatorestimatesparse
function/************************************************************************* This function estimates norm of the sparse M*N matrix A. INPUT PARAMETERS: State - norm estimator state, must be initialized with a call to NormEstimatorCreate() A - sparse M*N matrix, must be converted to CRS format prior to calling this function. After this function is over you can call NormEstimatorResults() to get estimate of the norm(A). -- ALGLIB -- Copyright 06.12.2011 by Bochkanov Sergey *************************************************************************/public static void alglib.normestimatorestimatesparse( normestimatorstate state, sparsematrix a)
normestimatorresults
function/************************************************************************* Matrix norm estimation results INPUT PARAMETERS: State - algorithm state OUTPUT PARAMETERS: Nrm - estimate of the matrix norm, Nrm>=0 -- ALGLIB -- Copyright 06.12.2011 by Bochkanov Sergey *************************************************************************/public static void alglib.normestimatorresults( normestimatorstate state, out double nrm)
normestimatorsetseed
function/************************************************************************* This function changes seed value used by algorithm. In some cases we need deterministic processing, i.e. subsequent calls must return equal results, in other cases we need non-deterministic algorithm which returns different results for the same matrix on every pass. Setting zero seed will lead to non-deterministic algorithm, while non-zero value will make our algorithm deterministic. INPUT PARAMETERS: State - norm estimator state, must be initialized with a call to NormEstimatorCreate() SeedVal - seed value, >=0. Zero value = non-deterministic algo. -- ALGLIB -- Copyright 06.12.2011 by Bochkanov Sergey *************************************************************************/public static void alglib.normestimatorsetseed( normestimatorstate state, int seedval)
odesolver
subpackage| odesolver_d1 | Solving y'=-y with ODE solver |
odesolverreport
class/************************************************************************* *************************************************************************/public class odesolverreport { public int nfev; public int terminationtype; }
odesolverstate
class/************************************************************************* *************************************************************************/public class odesolverstate { }
odesolverresults
function/************************************************************************* ODE solver results Called after OdeSolverIteration returned False. INPUT PARAMETERS: State - algorithm state (used by OdeSolverIteration). OUTPUT PARAMETERS: M - number of tabulated values, M>=1 XTbl - array[0..M-1], values of X YTbl - array[0..M-1,0..N-1], values of Y in X[i] Rep - solver report: * Rep.TerminationType completetion code: * -2 X is not ordered by ascending/descending or there are non-distinct X[], i.e. X[i]=X[i+1] * -1 incorrect parameters were specified * 1 task has been solved * Rep.NFEV contains number of function calculations -- ALGLIB -- Copyright 01.09.2009 by Bochkanov Sergey *************************************************************************/public static void alglib.odesolverresults( odesolverstate state, out int m, out double[] xtbl, out double[,] ytbl, out odesolverreport rep)
Examples: [1]
odesolverrkck
function/************************************************************************* Cash-Karp adaptive ODE solver. This subroutine solves ODE Y'=f(Y,x) with initial conditions Y(xs)=Ys (here Y may be single variable or vector of N variables). INPUT PARAMETERS: Y - initial conditions, array[0..N-1]. contains values of Y[] at X[0] N - system size X - points at which Y should be tabulated, array[0..M-1] integrations starts at X[0], ends at X[M-1], intermediate values at X[i] are returned too. SHOULD BE ORDERED BY ASCENDING OR BY DESCENDING!!!! M - number of intermediate points + first point + last point: * M>2 means that you need both Y(X[M-1]) and M-2 values at intermediate points * M=2 means that you want just to integrate from X[0] to X[1] and don't interested in intermediate values. * M=1 means that you don't want to integrate :) it is degenerate case, but it will be handled correctly. * M<1 means error Eps - tolerance (absolute/relative error on each step will be less than Eps). When passing: * Eps>0, it means desired ABSOLUTE error * Eps<0, it means desired RELATIVE error. Relative errors are calculated with respect to maximum values of Y seen so far. Be careful to use this criterion when starting from Y[] that are close to zero. H - initial step lenth, it will be adjusted automatically after the first step. If H=0, step will be selected automatically (usualy it will be equal to 0.001 of min(x[i]-x[j])). OUTPUT PARAMETERS State - structure which stores algorithm state between subsequent calls of OdeSolverIteration. Used for reverse communication. This structure should be passed to the OdeSolverIteration subroutine. SEE ALSO AutoGKSmoothW, AutoGKSingular, AutoGKIteration, AutoGKResults. -- ALGLIB -- Copyright 01.09.2009 by Bochkanov Sergey *************************************************************************/public static void alglib.odesolverrkck( double[] y, double[] x, double eps, double h, out odesolverstate state) public static void alglib.odesolverrkck( double[] y, int n, double[] x, int m, double eps, double h, out odesolverstate state)
Examples: [1]
odesolversolve
function/************************************************************************* This function is used to launcn iterations of ODE solver It accepts following parameters: diff - callback which calculates dy/dx for given y and x ptr - optional pointer which is passed to diff; can be NULL -- ALGLIB -- Copyright 01.09.2009 by Bochkanov Sergey *************************************************************************/public static void odesolversolve(odesolverstate state, ndimensional_ode_rp diff, object obj) { if( diff==null ) throw new alglibexception("ALGLIB: error in 'odesolversolve()' (diff is null)"); while( alglib.odesolveriteration(state) ) { if( state.needdy ) { diff(state.innerobj.y, state.innerobj.x, state.innerobj.dy, obj); continue; } throw new alglibexception("ALGLIB: unexpected error in 'odesolversolve'"); } }
Examples: [1]
public static void ode_function_1_diff(double[] y, double x, double[] dy, object obj) { // this callback calculates f(y[],x)=-y[0] dy[0] = -y[0]; } public static int Main(string[] args) { double[] y = new double[]{1}; double[] x = new double[]{0,1,2,3}; double eps = 0.00001; double h = 0; alglib.odesolverstate s; int m; double[] xtbl; double[,] ytbl; alglib.odesolverreport rep; alglib.odesolverrkck(y, x, eps, h, out s); alglib.odesolversolve(s, ode_function_1_diff, null); alglib.odesolverresults(s, out m, out xtbl, out ytbl, out rep); System.Console.WriteLine("{0}", m); // EXPECTED: 4 System.Console.WriteLine("{0}", alglib.ap.format(xtbl,2)); // EXPECTED: [0, 1, 2, 3] System.Console.WriteLine("{0}", alglib.ap.format(ytbl,2)); // EXPECTED: [[1], [0.367], [0.135], [0.050]] System.Console.ReadLine(); return 0; }
ortfac
subpackagecmatrixlq
function/************************************************************************* LQ decomposition of a rectangular complex matrix of size MxN Input parameters: A - matrix A whose indexes range within [0..M-1, 0..N-1] M - number of rows in matrix A. N - number of columns in matrix A. Output parameters: A - matrices Q and L in compact form Tau - array of scalar factors which are used to form matrix Q. Array whose indexes range within [0.. Min(M,N)-1] Matrix A is represented as A = LQ, where Q is an orthogonal matrix of size MxM, L - lower triangular (or lower trapezoid) matrix of size MxN. -- LAPACK routine (version 3.0) -- Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., Courant Institute, Argonne National Lab, and Rice University September 30, 1994 *************************************************************************/public static void alglib.cmatrixlq( ref complex[,] a, int m, int n, out complex[] tau)
cmatrixlqunpackl
function/************************************************************************* Unpacking of matrix L from the LQ decomposition of a matrix A Input parameters: A - matrices Q and L in compact form. Output of CMatrixLQ subroutine. M - number of rows in given matrix A. M>=0. N - number of columns in given matrix A. N>=0. Output parameters: L - matrix L, array[0..M-1, 0..N-1]. -- ALGLIB routine -- 17.02.2010 Bochkanov Sergey *************************************************************************/public static void alglib.cmatrixlqunpackl( complex[,] a, int m, int n, out complex[,] l)
cmatrixlqunpackq
function/************************************************************************* Partial unpacking of matrix Q from LQ decomposition of a complex matrix A. Input parameters: A - matrices Q and R in compact form. Output of CMatrixLQ subroutine . M - number of rows in matrix A. M>=0. N - number of columns in matrix A. N>=0. Tau - scalar factors which are used to form Q. Output of CMatrixLQ subroutine . QRows - required number of rows in matrix Q. N>=QColumns>=0. Output parameters: Q - first QRows rows of matrix Q. Array whose index ranges within [0..QRows-1, 0..N-1]. If QRows=0, array isn't changed. -- ALGLIB routine -- 17.02.2010 Bochkanov Sergey *************************************************************************/public static void alglib.cmatrixlqunpackq( complex[,] a, int m, int n, complex[] tau, int qrows, out complex[,] q)
cmatrixqr
function/************************************************************************* QR decomposition of a rectangular complex matrix of size MxN Input parameters: A - matrix A whose indexes range within [0..M-1, 0..N-1] M - number of rows in matrix A. N - number of columns in matrix A. Output parameters: A - matrices Q and R in compact form Tau - array of scalar factors which are used to form matrix Q. Array whose indexes range within [0.. Min(M,N)-1] Matrix A is represented as A = QR, where Q is an orthogonal matrix of size MxM, R - upper triangular (or upper trapezoid) matrix of size MxN. -- LAPACK routine (version 3.0) -- Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., Courant Institute, Argonne National Lab, and Rice University September 30, 1994 *************************************************************************/public static void alglib.cmatrixqr( ref complex[,] a, int m, int n, out complex[] tau)
cmatrixqrunpackq
function/************************************************************************* Partial unpacking of matrix Q from QR decomposition of a complex matrix A. Input parameters: A - matrices Q and R in compact form. Output of CMatrixQR subroutine . M - number of rows in matrix A. M>=0. N - number of columns in matrix A. N>=0. Tau - scalar factors which are used to form Q. Output of CMatrixQR subroutine . QColumns - required number of columns in matrix Q. M>=QColumns>=0. Output parameters: Q - first QColumns columns of matrix Q. Array whose index ranges within [0..M-1, 0..QColumns-1]. If QColumns=0, array isn't changed. -- ALGLIB routine -- 17.02.2010 Bochkanov Sergey *************************************************************************/public static void alglib.cmatrixqrunpackq( complex[,] a, int m, int n, complex[] tau, int qcolumns, out complex[,] q)
cmatrixqrunpackr
function/************************************************************************* Unpacking of matrix R from the QR decomposition of a matrix A Input parameters: A - matrices Q and R in compact form. Output of CMatrixQR subroutine. M - number of rows in given matrix A. M>=0. N - number of columns in given matrix A. N>=0. Output parameters: R - matrix R, array[0..M-1, 0..N-1]. -- ALGLIB routine -- 17.02.2010 Bochkanov Sergey *************************************************************************/public static void alglib.cmatrixqrunpackr( complex[,] a, int m, int n, out complex[,] r)
hmatrixtd
function/************************************************************************* Reduction of a Hermitian matrix which is given by its higher or lower triangular part to a real tridiagonal matrix using unitary similarity transformation: Q'*A*Q = T. Input parameters: A - matrix to be transformed array with elements [0..N-1, 0..N-1]. N - size of matrix A. IsUpper - storage format. If IsUpper = True, then matrix A is given by its upper triangle, and the lower triangle is not used and not modified by the algorithm, and vice versa if IsUpper = False. Output parameters: A - matrices T and Q in compact form (see lower) Tau - array of factors which are forming matrices H(i) array with elements [0..N-2]. D - main diagonal of real symmetric matrix T. array with elements [0..N-1]. E - secondary diagonal of real symmetric matrix T. array with elements [0..N-2]. If IsUpper=True, the matrix Q is represented as a product of elementary reflectors Q = H(n-2) . . . H(2) H(0). Each H(i) has the form H(i) = I - tau * v * v' where tau is a complex scalar, and v is a complex vector with v(i+1:n-1) = 0, v(i) = 1, v(0:i-1) is stored on exit in A(0:i-1,i+1), and tau in TAU(i). If IsUpper=False, the matrix Q is represented as a product of elementary reflectors Q = H(0) H(2) . . . H(n-2). Each H(i) has the form H(i) = I - tau * v * v' where tau is a complex scalar, and v is a complex vector with v(0:i) = 0, v(i+1) = 1, v(i+2:n-1) is stored on exit in A(i+2:n-1,i), and tau in TAU(i). The contents of A on exit are illustrated by the following examples with n = 5: if UPLO = 'U': if UPLO = 'L': ( d e v1 v2 v3 ) ( d ) ( d e v2 v3 ) ( e d ) ( d e v3 ) ( v0 e d ) ( d e ) ( v0 v1 e d ) ( d ) ( v0 v1 v2 e d ) where d and e denote diagonal and off-diagonal elements of T, and vi denotes an element of the vector defining H(i). -- LAPACK routine (version 3.0) -- Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., Courant Institute, Argonne National Lab, and Rice University October 31, 1992 *************************************************************************/public static void alglib.hmatrixtd( ref complex[,] a, int n, bool isupper, out complex[] tau, out double[] d, out double[] e)
hmatrixtdunpackq
function/************************************************************************* Unpacking matrix Q which reduces a Hermitian matrix to a real tridiagonal form. Input parameters: A - the result of a HMatrixTD subroutine N - size of matrix A. IsUpper - storage format (a parameter of HMatrixTD subroutine) Tau - the result of a HMatrixTD subroutine Output parameters: Q - transformation matrix. array with elements [0..N-1, 0..N-1]. -- ALGLIB -- Copyright 2005-2010 by Bochkanov Sergey *************************************************************************/public static void alglib.hmatrixtdunpackq( complex[,] a, int n, bool isupper, complex[] tau, out complex[,] q)
rmatrixbd
function/************************************************************************* Reduction of a rectangular matrix to bidiagonal form The algorithm reduces the rectangular matrix A to bidiagonal form by orthogonal transformations P and Q: A = Q*B*P. Input parameters: A - source matrix. array[0..M-1, 0..N-1] M - number of rows in matrix A. N - number of columns in matrix A. Output parameters: A - matrices Q, B, P in compact form (see below). TauQ - scalar factors which are used to form matrix Q. TauP - scalar factors which are used to form matrix P. The main diagonal and one of the secondary diagonals of matrix A are replaced with bidiagonal matrix B. Other elements contain elementary reflections which form MxM matrix Q and NxN matrix P, respectively. If M>=N, B is the upper bidiagonal MxN matrix and is stored in the corresponding elements of matrix A. Matrix Q is represented as a product of elementary reflections Q = H(0)*H(1)*...*H(n-1), where H(i) = 1-tau*v*v'. Here tau is a scalar which is stored in TauQ[i], and vector v has the following structure: v(0:i-1)=0, v(i)=1, v(i+1:m-1) is stored in elements A(i+1:m-1,i). Matrix P is as follows: P = G(0)*G(1)*...*G(n-2), where G(i) = 1 - tau*u*u'. Tau is stored in TauP[i], u(0:i)=0, u(i+1)=1, u(i+2:n-1) is stored in elements A(i,i+2:n-1). If M<N, B is the lower bidiagonal MxN matrix and is stored in the corresponding elements of matrix A. Q = H(0)*H(1)*...*H(m-2), where H(i) = 1 - tau*v*v', tau is stored in TauQ, v(0:i)=0, v(i+1)=1, v(i+2:m-1) is stored in elements A(i+2:m-1,i). P = G(0)*G(1)*...*G(m-1), G(i) = 1-tau*u*u', tau is stored in TauP, u(0:i-1)=0, u(i)=1, u(i+1:n-1) is stored in A(i,i+1:n-1). EXAMPLE: m=6, n=5 (m > n): m=5, n=6 (m < n): ( d e u1 u1 u1 ) ( d u1 u1 u1 u1 u1 ) ( v1 d e u2 u2 ) ( e d u2 u2 u2 u2 ) ( v1 v2 d e u3 ) ( v1 e d u3 u3 u3 ) ( v1 v2 v3 d e ) ( v1 v2 e d u4 u4 ) ( v1 v2 v3 v4 d ) ( v1 v2 v3 e d u5 ) ( v1 v2 v3 v4 v5 ) Here vi and ui are vectors which form H(i) and G(i), and d and e - are the diagonal and off-diagonal elements of matrix B. -- LAPACK routine (version 3.0) -- Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., Courant Institute, Argonne National Lab, and Rice University September 30, 1994. Sergey Bochkanov, ALGLIB project, translation from FORTRAN to pseudocode, 2007-2010. *************************************************************************/public static void alglib.rmatrixbd( ref double[,] a, int m, int n, out double[] tauq, out double[] taup)
rmatrixbdmultiplybyp
function/************************************************************************* Multiplication by matrix P which reduces matrix A to bidiagonal form. The algorithm allows pre- or post-multiply by P or P'. Input parameters: QP - matrices Q and P in compact form. Output of RMatrixBD subroutine. M - number of rows in matrix A. N - number of columns in matrix A. TAUP - scalar factors which are used to form P. Output of RMatrixBD subroutine. Z - multiplied matrix. Array whose indexes range within [0..ZRows-1,0..ZColumns-1]. ZRows - number of rows in matrix Z. If FromTheRight=False, ZRows=N, otherwise ZRows can be arbitrary. ZColumns - number of columns in matrix Z. If FromTheRight=True, ZColumns=N, otherwise ZColumns can be arbitrary. FromTheRight - pre- or post-multiply. DoTranspose - multiply by P or P'. Output parameters: Z - product of Z and P. Array whose indexes range within [0..ZRows-1,0..ZColumns-1]. If ZRows=0 or ZColumns=0, the array is not modified. -- ALGLIB -- 2005-2010 Bochkanov Sergey *************************************************************************/public static void alglib.rmatrixbdmultiplybyp( double[,] qp, int m, int n, double[] taup, ref double[,] z, int zrows, int zcolumns, bool fromtheright, bool dotranspose)
rmatrixbdmultiplybyq
function/************************************************************************* Multiplication by matrix Q which reduces matrix A to bidiagonal form. The algorithm allows pre- or post-multiply by Q or Q'. Input parameters: QP - matrices Q and P in compact form. Output of ToBidiagonal subroutine. M - number of rows in matrix A. N - number of columns in matrix A. TAUQ - scalar factors which are used to form Q. Output of ToBidiagonal subroutine. Z - multiplied matrix. array[0..ZRows-1,0..ZColumns-1] ZRows - number of rows in matrix Z. If FromTheRight=False, ZRows=M, otherwise ZRows can be arbitrary. ZColumns - number of columns in matrix Z. If FromTheRight=True, ZColumns=M, otherwise ZColumns can be arbitrary. FromTheRight - pre- or post-multiply. DoTranspose - multiply by Q or Q'. Output parameters: Z - product of Z and Q. Array[0..ZRows-1,0..ZColumns-1] If ZRows=0 or ZColumns=0, the array is not modified. -- ALGLIB -- 2005-2010 Bochkanov Sergey *************************************************************************/public static void alglib.rmatrixbdmultiplybyq( double[,] qp, int m, int n, double[] tauq, ref double[,] z, int zrows, int zcolumns, bool fromtheright, bool dotranspose)
rmatrixbdunpackdiagonals
function/************************************************************************* Unpacking of the main and secondary diagonals of bidiagonal decomposition of matrix A. Input parameters: B - output of RMatrixBD subroutine. M - number of rows in matrix B. N - number of columns in matrix B. Output parameters: IsUpper - True, if the matrix is upper bidiagonal. otherwise IsUpper is False. D - the main diagonal. Array whose index ranges within [0..Min(M,N)-1]. E - the secondary diagonal (upper or lower, depending on the value of IsUpper). Array index ranges within [0..Min(M,N)-1], the last element is not used. -- ALGLIB -- 2005-2010 Bochkanov Sergey *************************************************************************/public static void alglib.rmatrixbdunpackdiagonals( double[,] b, int m, int n, out bool isupper, out double[] d, out double[] e)
rmatrixbdunpackpt
function/************************************************************************* Unpacking matrix P which reduces matrix A to bidiagonal form. The subroutine returns transposed matrix P. Input parameters: QP - matrices Q and P in compact form. Output of ToBidiagonal subroutine. M - number of rows in matrix A. N - number of columns in matrix A. TAUP - scalar factors which are used to form P. Output of ToBidiagonal subroutine. PTRows - required number of rows of matrix P^T. N >= PTRows >= 0. Output parameters: PT - first PTRows columns of matrix P^T Array[0..PTRows-1, 0..N-1] If PTRows=0, the array is not modified. -- ALGLIB -- 2005-2010 Bochkanov Sergey *************************************************************************/public static void alglib.rmatrixbdunpackpt( double[,] qp, int m, int n, double[] taup, int ptrows, out double[,] pt)
rmatrixbdunpackq
function/************************************************************************* Unpacking matrix Q which reduces a matrix to bidiagonal form. Input parameters: QP - matrices Q and P in compact form. Output of ToBidiagonal subroutine. M - number of rows in matrix A. N - number of columns in matrix A. TAUQ - scalar factors which are used to form Q. Output of ToBidiagonal subroutine. QColumns - required number of columns in matrix Q. M>=QColumns>=0. Output parameters: Q - first QColumns columns of matrix Q. Array[0..M-1, 0..QColumns-1] If QColumns=0, the array is not modified. -- ALGLIB -- 2005-2010 Bochkanov Sergey *************************************************************************/public static void alglib.rmatrixbdunpackq( double[,] qp, int m, int n, double[] tauq, int qcolumns, out double[,] q)
rmatrixhessenberg
function/************************************************************************* Reduction of a square matrix to upper Hessenberg form: Q'*A*Q = H, where Q is an orthogonal matrix, H - Hessenberg matrix. Input parameters: A - matrix A with elements [0..N-1, 0..N-1] N - size of matrix A. Output parameters: A - matrices Q and P in compact form (see below). Tau - array of scalar factors which are used to form matrix Q. Array whose index ranges within [0..N-2] Matrix H is located on the main diagonal, on the lower secondary diagonal and above the main diagonal of matrix A. The elements which are used to form matrix Q are situated in array Tau and below the lower secondary diagonal of matrix A as follows: Matrix Q is represented as a product of elementary reflections Q = H(0)*H(2)*...*H(n-2), where each H(i) is given by H(i) = 1 - tau * v * (v^T) where tau is a scalar stored in Tau[I]; v - is a real vector, so that v(0:i) = 0, v(i+1) = 1, v(i+2:n-1) stored in A(i+2:n-1,i). -- LAPACK routine (version 3.0) -- Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., Courant Institute, Argonne National Lab, and Rice University October 31, 1992 *************************************************************************/public static void alglib.rmatrixhessenberg( ref double[,] a, int n, out double[] tau)
rmatrixhessenbergunpackh
function/************************************************************************* Unpacking matrix H (the result of matrix A reduction to upper Hessenberg form) Input parameters: A - output of RMatrixHessenberg subroutine. N - size of matrix A. Output parameters: H - matrix H. Array whose indexes range within [0..N-1, 0..N-1]. -- ALGLIB -- 2005-2010 Bochkanov Sergey *************************************************************************/public static void alglib.rmatrixhessenbergunpackh( double[,] a, int n, out double[,] h)
rmatrixhessenbergunpackq
function/************************************************************************* Unpacking matrix Q which reduces matrix A to upper Hessenberg form Input parameters: A - output of RMatrixHessenberg subroutine. N - size of matrix A. Tau - scalar factors which are used to form Q. Output of RMatrixHessenberg subroutine. Output parameters: Q - matrix Q. Array whose indexes range within [0..N-1, 0..N-1]. -- ALGLIB -- 2005-2010 Bochkanov Sergey *************************************************************************/public static void alglib.rmatrixhessenbergunpackq( double[,] a, int n, double[] tau, out double[,] q)
rmatrixlq
function/************************************************************************* LQ decomposition of a rectangular matrix of size MxN Input parameters: A - matrix A whose indexes range within [0..M-1, 0..N-1]. M - number of rows in matrix A. N - number of columns in matrix A. Output parameters: A - matrices L and Q in compact form (see below) Tau - array of scalar factors which are used to form matrix Q. Array whose index ranges within [0..Min(M,N)-1]. Matrix A is represented as A = LQ, where Q is an orthogonal matrix of size MxM, L - lower triangular (or lower trapezoid) matrix of size M x N. The elements of matrix L are located on and below the main diagonal of matrix A. The elements which are located in Tau array and above the main diagonal of matrix A are used to form matrix Q as follows: Matrix Q is represented as a product of elementary reflections Q = H(k-1)*H(k-2)*...*H(1)*H(0), where k = min(m,n), and each H(i) is of the form H(i) = 1 - tau * v * (v^T) where tau is a scalar stored in Tau[I]; v - real vector, so that v(0:i-1)=0, v(i) = 1, v(i+1:n-1) stored in A(i,i+1:n-1). -- ALGLIB routine -- 17.02.2010 Bochkanov Sergey *************************************************************************/public static void alglib.rmatrixlq( ref double[,] a, int m, int n, out double[] tau)
rmatrixlqunpackl
function/************************************************************************* Unpacking of matrix L from the LQ decomposition of a matrix A Input parameters: A - matrices Q and L in compact form. Output of RMatrixLQ subroutine. M - number of rows in given matrix A. M>=0. N - number of columns in given matrix A. N>=0. Output parameters: L - matrix L, array[0..M-1, 0..N-1]. -- ALGLIB routine -- 17.02.2010 Bochkanov Sergey *************************************************************************/public static void alglib.rmatrixlqunpackl( double[,] a, int m, int n, out double[,] l)
rmatrixlqunpackq
function/************************************************************************* Partial unpacking of matrix Q from the LQ decomposition of a matrix A Input parameters: A - matrices L and Q in compact form. Output of RMatrixLQ subroutine. M - number of rows in given matrix A. M>=0. N - number of columns in given matrix A. N>=0. Tau - scalar factors which are used to form Q. Output of the RMatrixLQ subroutine. QRows - required number of rows in matrix Q. N>=QRows>=0. Output parameters: Q - first QRows rows of matrix Q. Array whose indexes range within [0..QRows-1, 0..N-1]. If QRows=0, the array remains unchanged. -- ALGLIB routine -- 17.02.2010 Bochkanov Sergey *************************************************************************/public static void alglib.rmatrixlqunpackq( double[,] a, int m, int n, double[] tau, int qrows, out double[,] q)
rmatrixqr
function/************************************************************************* QR decomposition of a rectangular matrix of size MxN Input parameters: A - matrix A whose indexes range within [0..M-1, 0..N-1]. M - number of rows in matrix A. N - number of columns in matrix A. Output parameters: A - matrices Q and R in compact form (see below). Tau - array of scalar factors which are used to form matrix Q. Array whose index ranges within [0.. Min(M-1,N-1)]. Matrix A is represented as A = QR, where Q is an orthogonal matrix of size MxM, R - upper triangular (or upper trapezoid) matrix of size M x N. The elements of matrix R are located on and above the main diagonal of matrix A. The elements which are located in Tau array and below the main diagonal of matrix A are used to form matrix Q as follows: Matrix Q is represented as a product of elementary reflections Q = H(0)*H(2)*...*H(k-1), where k = min(m,n), and each H(i) is in the form H(i) = 1 - tau * v * (v^T) where tau is a scalar stored in Tau[I]; v - real vector, so that v(0:i-1) = 0, v(i) = 1, v(i+1:m-1) stored in A(i+1:m-1,i). -- ALGLIB routine -- 17.02.2010 Bochkanov Sergey *************************************************************************/public static void alglib.rmatrixqr( ref double[,] a, int m, int n, out double[] tau)
rmatrixqrunpackq
function/************************************************************************* Partial unpacking of matrix Q from the QR decomposition of a matrix A Input parameters: A - matrices Q and R in compact form. Output of RMatrixQR subroutine. M - number of rows in given matrix A. M>=0. N - number of columns in given matrix A. N>=0. Tau - scalar factors which are used to form Q. Output of the RMatrixQR subroutine. QColumns - required number of columns of matrix Q. M>=QColumns>=0. Output parameters: Q - first QColumns columns of matrix Q. Array whose indexes range within [0..M-1, 0..QColumns-1]. If QColumns=0, the array remains unchanged. -- ALGLIB routine -- 17.02.2010 Bochkanov Sergey *************************************************************************/public static void alglib.rmatrixqrunpackq( double[,] a, int m, int n, double[] tau, int qcolumns, out double[,] q)
rmatrixqrunpackr
function/************************************************************************* Unpacking of matrix R from the QR decomposition of a matrix A Input parameters: A - matrices Q and R in compact form. Output of RMatrixQR subroutine. M - number of rows in given matrix A. M>=0. N - number of columns in given matrix A. N>=0. Output parameters: R - matrix R, array[0..M-1, 0..N-1]. -- ALGLIB routine -- 17.02.2010 Bochkanov Sergey *************************************************************************/public static void alglib.rmatrixqrunpackr( double[,] a, int m, int n, out double[,] r)
smatrixtd
function/************************************************************************* Reduction of a symmetric matrix which is given by its higher or lower triangular part to a tridiagonal matrix using orthogonal similarity transformation: Q'*A*Q=T. Input parameters: A - matrix to be transformed array with elements [0..N-1, 0..N-1]. N - size of matrix A. IsUpper - storage format. If IsUpper = True, then matrix A is given by its upper triangle, and the lower triangle is not used and not modified by the algorithm, and vice versa if IsUpper = False. Output parameters: A - matrices T and Q in compact form (see lower) Tau - array of factors which are forming matrices H(i) array with elements [0..N-2]. D - main diagonal of symmetric matrix T. array with elements [0..N-1]. E - secondary diagonal of symmetric matrix T. array with elements [0..N-2]. If IsUpper=True, the matrix Q is represented as a product of elementary reflectors Q = H(n-2) . . . H(2) H(0). Each H(i) has the form H(i) = I - tau * v * v' where tau is a real scalar, and v is a real vector with v(i+1:n-1) = 0, v(i) = 1, v(0:i-1) is stored on exit in A(0:i-1,i+1), and tau in TAU(i). If IsUpper=False, the matrix Q is represented as a product of elementary reflectors Q = H(0) H(2) . . . H(n-2). Each H(i) has the form H(i) = I - tau * v * v' where tau is a real scalar, and v is a real vector with v(0:i) = 0, v(i+1) = 1, v(i+2:n-1) is stored on exit in A(i+2:n-1,i), and tau in TAU(i). The contents of A on exit are illustrated by the following examples with n = 5: if UPLO = 'U': if UPLO = 'L': ( d e v1 v2 v3 ) ( d ) ( d e v2 v3 ) ( e d ) ( d e v3 ) ( v0 e d ) ( d e ) ( v0 v1 e d ) ( d ) ( v0 v1 v2 e d ) where d and e denote diagonal and off-diagonal elements of T, and vi denotes an element of the vector defining H(i). -- LAPACK routine (version 3.0) -- Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., Courant Institute, Argonne National Lab, and Rice University October 31, 1992 *************************************************************************/public static void alglib.smatrixtd( ref double[,] a, int n, bool isupper, out double[] tau, out double[] d, out double[] e)
smatrixtdunpackq
function/************************************************************************* Unpacking matrix Q which reduces symmetric matrix to a tridiagonal form. Input parameters: A - the result of a SMatrixTD subroutine N - size of matrix A. IsUpper - storage format (a parameter of SMatrixTD subroutine) Tau - the result of a SMatrixTD subroutine Output parameters: Q - transformation matrix. array with elements [0..N-1, 0..N-1]. -- ALGLIB -- Copyright 2005-2010 by Bochkanov Sergey *************************************************************************/public static void alglib.smatrixtdunpackq( double[,] a, int n, bool isupper, double[] tau, out double[,] q)
pca
subpackagepcabuildbasis
function/************************************************************************* Principal components analysis Subroutine builds orthogonal basis where first axis corresponds to direction with maximum variance, second axis maximizes variance in subspace orthogonal to first axis and so on. It should be noted that, unlike LDA, PCA does not use class labels. INPUT PARAMETERS: X - dataset, array[0..NPoints-1,0..NVars-1]. matrix contains ONLY INDEPENDENT VARIABLES. NPoints - dataset size, NPoints>=0 NVars - number of independent variables, NVars>=1 OUTPUT PARAMETERS: Info - return code: * -4, if SVD subroutine haven't converged * -1, if wrong parameters has been passed (NPoints<0, NVars<1) * 1, if task is solved S2 - array[0..NVars-1]. variance values corresponding to basis vectors. V - array[0..NVars-1,0..NVars-1] matrix, whose columns store basis vectors. -- ALGLIB -- Copyright 25.08.2008 by Bochkanov Sergey *************************************************************************/public static void alglib.pcabuildbasis( double[,] x, int npoints, int nvars, out int info, out double[] s2, out double[,] v)
poissondistr
subpackageinvpoissondistribution
function/************************************************************************* Inverse Poisson distribution Finds the Poisson variable x such that the integral from 0 to x of the Poisson density is equal to the given probability y. This is accomplished using the inverse gamma integral function and the relation m = igami( k+1, y ). ACCURACY: See inverse incomplete gamma function Cephes Math Library Release 2.8: June, 2000 Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier *************************************************************************/public static double alglib.invpoissondistribution(int k, double y)
poissoncdistribution
function/************************************************************************* Complemented Poisson distribution Returns the sum of the terms k+1 to infinity of the Poisson distribution: inf. j -- -m m > e -- -- j! j=k+1 The terms are not summed directly; instead the incomplete gamma integral is employed, according to the formula y = pdtrc( k, m ) = igam( k+1, m ). The arguments must both be positive. ACCURACY: See incomplete gamma function Cephes Math Library Release 2.8: June, 2000 Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier *************************************************************************/public static double alglib.poissoncdistribution(int k, double m)
poissondistribution
function/************************************************************************* Poisson distribution Returns the sum of the first k+1 terms of the Poisson distribution: k j -- -m m > e -- -- j! j=0 The terms are not summed directly; instead the incomplete gamma integral is employed, according to the relation y = pdtr( k, m ) = igamc( k+1, m ). The arguments must both be positive. ACCURACY: See incomplete gamma function Cephes Math Library Release 2.8: June, 2000 Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier *************************************************************************/public static double alglib.poissondistribution(int k, double m)
polint
subpackage| polint_d_calcdiff | Interpolation and differentiation using barycentric representation | |
| polint_d_conv | Conversion between power basis and barycentric representation | |
| polint_d_spec | Polynomial interpolation on special grids (equidistant, Chebyshev I/II) |
polynomialbar2cheb
function/************************************************************************* Conversion from barycentric representation to Chebyshev basis. This function has O(N^2) complexity. INPUT PARAMETERS: P - polynomial in barycentric form A,B - base interval for Chebyshev polynomials (see below) A<>B OUTPUT PARAMETERS T - coefficients of Chebyshev representation; P(x) = sum { T[i]*Ti(2*(x-A)/(B-A)-1), i=0..N-1 }, where Ti - I-th Chebyshev polynomial. NOTES: barycentric interpolant passed as P may be either polynomial obtained from polynomial interpolation/ fitting or rational function which is NOT polynomial. We can't distinguish between these two cases, and this algorithm just tries to work assuming that P IS a polynomial. If not, algorithm will return results, but they won't have any meaning. -- ALGLIB -- Copyright 30.09.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.polynomialbar2cheb( barycentricinterpolant p, double a, double b, out double[] t)
polynomialbar2pow
function/************************************************************************* Conversion from barycentric representation to power basis. This function has O(N^2) complexity. INPUT PARAMETERS: P - polynomial in barycentric form C - offset (see below); 0.0 is used as default value. S - scale (see below); 1.0 is used as default value. S<>0. OUTPUT PARAMETERS A - coefficients, P(x) = sum { A[i]*((X-C)/S)^i, i=0..N-1 } N - number of coefficients (polynomial degree plus 1) NOTES: 1. this function accepts offset and scale, which can be set to improve numerical properties of polynomial. For example, if P was obtained as result of interpolation on [-1,+1], you can set C=0 and S=1 and represent P as sum of 1, x, x^2, x^3 and so on. In most cases you it is exactly what you need. However, if your interpolation model was built on [999,1001], you will see significant growth of numerical errors when using {1, x, x^2, x^3} as basis. Representing P as sum of 1, (x-1000), (x-1000)^2, (x-1000)^3 will be better option. Such representation can be obtained by using 1000.0 as offset C and 1.0 as scale S. 2. power basis is ill-conditioned and tricks described above can't solve this problem completely. This function will return coefficients in any case, but for N>8 they will become unreliable. However, N's less than 5 are pretty safe. 3. barycentric interpolant passed as P may be either polynomial obtained from polynomial interpolation/ fitting or rational function which is NOT polynomial. We can't distinguish between these two cases, and this algorithm just tries to work assuming that P IS a polynomial. If not, algorithm will return results, but they won't have any meaning. -- ALGLIB -- Copyright 30.09.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.polynomialbar2pow( barycentricinterpolant p, out double[] a) public static void alglib.polynomialbar2pow( barycentricinterpolant p, double c, double s, out double[] a)
Examples: [1]
polynomialbuild
function/************************************************************************* Lagrange intepolant: generation of the model on the general grid. This function has O(N^2) complexity. INPUT PARAMETERS: X - abscissas, array[0..N-1] Y - function values, array[0..N-1] N - number of points, N>=1 OUTPUT PARAMETERS P - barycentric model which represents Lagrange interpolant (see ratint unit info and BarycentricCalc() description for more information). -- ALGLIB -- Copyright 02.12.2009 by Bochkanov Sergey *************************************************************************/public static void alglib.polynomialbuild( double[] x, double[] y, out barycentricinterpolant p) public static void alglib.polynomialbuild( double[] x, double[] y, int n, out barycentricinterpolant p)
Examples: [1]
polynomialbuildcheb1
function/************************************************************************* Lagrange intepolant on Chebyshev grid (first kind). This function has O(N) complexity. INPUT PARAMETERS: A - left boundary of [A,B] B - right boundary of [A,B] Y - function values at the nodes, array[0..N-1], Y[I] = Y(0.5*(B+A) + 0.5*(B-A)*Cos(PI*(2*i+1)/(2*n))) N - number of points, N>=1 for N=1 a constant model is constructed. OUTPUT PARAMETERS P - barycentric model which represents Lagrange interpolant (see ratint unit info and BarycentricCalc() description for more information). -- ALGLIB -- Copyright 03.12.2009 by Bochkanov Sergey *************************************************************************/public static void alglib.polynomialbuildcheb1( double a, double b, double[] y, out barycentricinterpolant p) public static void alglib.polynomialbuildcheb1( double a, double b, double[] y, int n, out barycentricinterpolant p)
Examples: [1]
polynomialbuildcheb2
function/************************************************************************* Lagrange intepolant on Chebyshev grid (second kind). This function has O(N) complexity. INPUT PARAMETERS: A - left boundary of [A,B] B - right boundary of [A,B] Y - function values at the nodes, array[0..N-1], Y[I] = Y(0.5*(B+A) + 0.5*(B-A)*Cos(PI*i/(n-1))) N - number of points, N>=1 for N=1 a constant model is constructed. OUTPUT PARAMETERS P - barycentric model which represents Lagrange interpolant (see ratint unit info and BarycentricCalc() description for more information). -- ALGLIB -- Copyright 03.12.2009 by Bochkanov Sergey *************************************************************************/public static void alglib.polynomialbuildcheb2( double a, double b, double[] y, out barycentricinterpolant p) public static void alglib.polynomialbuildcheb2( double a, double b, double[] y, int n, out barycentricinterpolant p)
Examples: [1]
polynomialbuildeqdist
function/************************************************************************* Lagrange intepolant: generation of the model on equidistant grid. This function has O(N) complexity. INPUT PARAMETERS: A - left boundary of [A,B] B - right boundary of [A,B] Y - function values at the nodes, array[0..N-1] N - number of points, N>=1 for N=1 a constant model is constructed. OUTPUT PARAMETERS P - barycentric model which represents Lagrange interpolant (see ratint unit info and BarycentricCalc() description for more information). -- ALGLIB -- Copyright 03.12.2009 by Bochkanov Sergey *************************************************************************/public static void alglib.polynomialbuildeqdist( double a, double b, double[] y, out barycentricinterpolant p) public static void alglib.polynomialbuildeqdist( double a, double b, double[] y, int n, out barycentricinterpolant p)
Examples: [1]
polynomialcalccheb1
function/************************************************************************* Fast polynomial interpolation function on Chebyshev points (first kind) with O(N) complexity. INPUT PARAMETERS: A - left boundary of [A,B] B - right boundary of [A,B] F - function values, array[0..N-1] N - number of points on Chebyshev grid (first kind), X[i] = 0.5*(B+A) + 0.5*(B-A)*Cos(PI*(2*i+1)/(2*n)) for N=1 a constant model is constructed. T - position where P(x) is calculated RESULT value of the Lagrange interpolant at T IMPORTANT this function provides fast interface which is not overflow-safe nor it is very precise. the best option is to use PolIntBuildCheb1()/BarycentricCalc() subroutines unless you are pretty sure that your data will not result in overflow. -- ALGLIB -- Copyright 02.12.2009 by Bochkanov Sergey *************************************************************************/public static double alglib.polynomialcalccheb1( double a, double b, double[] f, double t) public static double alglib.polynomialcalccheb1( double a, double b, double[] f, int n, double t)
Examples: [1]
polynomialcalccheb2
function/************************************************************************* Fast polynomial interpolation function on Chebyshev points (second kind) with O(N) complexity. INPUT PARAMETERS: A - left boundary of [A,B] B - right boundary of [A,B] F - function values, array[0..N-1] N - number of points on Chebyshev grid (second kind), X[i] = 0.5*(B+A) + 0.5*(B-A)*Cos(PI*i/(n-1)) for N=1 a constant model is constructed. T - position where P(x) is calculated RESULT value of the Lagrange interpolant at T IMPORTANT this function provides fast interface which is not overflow-safe nor it is very precise. the best option is to use PolIntBuildCheb2()/BarycentricCalc() subroutines unless you are pretty sure that your data will not result in overflow. -- ALGLIB -- Copyright 02.12.2009 by Bochkanov Sergey *************************************************************************/public static double alglib.polynomialcalccheb2( double a, double b, double[] f, double t) public static double alglib.polynomialcalccheb2( double a, double b, double[] f, int n, double t)
Examples: [1]
polynomialcalceqdist
function/************************************************************************* Fast equidistant polynomial interpolation function with O(N) complexity INPUT PARAMETERS: A - left boundary of [A,B] B - right boundary of [A,B] F - function values, array[0..N-1] N - number of points on equidistant grid, N>=1 for N=1 a constant model is constructed. T - position where P(x) is calculated RESULT value of the Lagrange interpolant at T IMPORTANT this function provides fast interface which is not overflow-safe nor it is very precise. the best option is to use PolynomialBuildEqDist()/BarycentricCalc() subroutines unless you are pretty sure that your data will not result in overflow. -- ALGLIB -- Copyright 02.12.2009 by Bochkanov Sergey *************************************************************************/public static double alglib.polynomialcalceqdist( double a, double b, double[] f, double t) public static double alglib.polynomialcalceqdist( double a, double b, double[] f, int n, double t)
Examples: [1]
polynomialcheb2bar
function/************************************************************************* Conversion from Chebyshev basis to barycentric representation. This function has O(N^2) complexity. INPUT PARAMETERS: T - coefficients of Chebyshev representation; P(x) = sum { T[i]*Ti(2*(x-A)/(B-A)-1), i=0..N }, where Ti - I-th Chebyshev polynomial. N - number of coefficients: * if given, only leading N elements of T are used * if not given, automatically determined from size of T A,B - base interval for Chebyshev polynomials (see above) A<B OUTPUT PARAMETERS P - polynomial in barycentric form -- ALGLIB -- Copyright 30.09.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.polynomialcheb2bar( double[] t, double a, double b, out barycentricinterpolant p) public static void alglib.polynomialcheb2bar( double[] t, int n, double a, double b, out barycentricinterpolant p)
polynomialpow2bar
function/************************************************************************* Conversion from power basis to barycentric representation. This function has O(N^2) complexity. INPUT PARAMETERS: A - coefficients, P(x) = sum { A[i]*((X-C)/S)^i, i=0..N-1 } N - number of coefficients (polynomial degree plus 1) * if given, only leading N elements of A are used * if not given, automatically determined from size of A C - offset (see below); 0.0 is used as default value. S - scale (see below); 1.0 is used as default value. S<>0. OUTPUT PARAMETERS P - polynomial in barycentric form NOTES: 1. this function accepts offset and scale, which can be set to improve numerical properties of polynomial. For example, if you interpolate on [-1,+1], you can set C=0 and S=1 and convert from sum of 1, x, x^2, x^3 and so on. In most cases you it is exactly what you need. However, if your interpolation model was built on [999,1001], you will see significant growth of numerical errors when using {1, x, x^2, x^3} as input basis. Converting from sum of 1, (x-1000), (x-1000)^2, (x-1000)^3 will be better option (you have to specify 1000.0 as offset C and 1.0 as scale S). 2. power basis is ill-conditioned and tricks described above can't solve this problem completely. This function will return barycentric model in any case, but for N>8 accuracy well degrade. However, N's less than 5 are pretty safe. -- ALGLIB -- Copyright 30.09.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.polynomialpow2bar( double[] a, out barycentricinterpolant p) public static void alglib.polynomialpow2bar( double[] a, int n, double c, double s, out barycentricinterpolant p)
Examples: [1]
public static int Main(string[] args) { // // Here we demonstrate polynomial interpolation and differentiation // of y=x^2-x sampled at [0,1,2]. Barycentric representation of polynomial is used. // double[] x = new double[]{0,1,2}; double[] y = new double[]{0,0,2}; double t = -1; double v; double dv; double d2v; alglib.barycentricinterpolant p; // barycentric model is created alglib.polynomialbuild(x, y, out p); // barycentric interpolation is demonstrated v = alglib.barycentriccalc(p, t); System.Console.WriteLine("{0:F4", v); // EXPECTED: 2.0 // barycentric differentation is demonstrated alglib.barycentricdiff1(p, t, out v, out dv); System.Console.WriteLine("{0:F4", v); // EXPECTED: 2.0 System.Console.WriteLine("{0:F4", dv); // EXPECTED: -3.0 // second derivatives with barycentric representation alglib.barycentricdiff1(p, t, out v, out dv); System.Console.WriteLine("{0:F4", v); // EXPECTED: 2.0 System.Console.WriteLine("{0:F4", dv); // EXPECTED: -3.0 alglib.barycentricdiff2(p, t, out v, out dv, out d2v); System.Console.WriteLine("{0:F4", v); // EXPECTED: 2.0 System.Console.WriteLine("{0:F4", dv); // EXPECTED: -3.0 System.Console.WriteLine("{0:F4", d2v); // EXPECTED: 2.0 System.Console.ReadLine(); return 0; }
public static int Main(string[] args) { // // Here we demonstrate conversion of y=x^2-x // between power basis and barycentric representation. // double[] a = new double[]{0,-1,+1}; double t = 2; double[] a2; double v; alglib.barycentricinterpolant p; // // a=[0,-1,+1] is decomposition of y=x^2-x in the power basis: // // y = 0 - 1*x + 1*x^2 // // We convert it to the barycentric form. // alglib.polynomialpow2bar(a, out p); // now we have barycentric interpolation; we can use it for interpolation v = alglib.barycentriccalc(p, t); System.Console.WriteLine("{0:F2", v); // EXPECTED: 2.0 // we can also convert back from barycentric representation to power basis alglib.polynomialbar2pow(p, out a2); System.Console.WriteLine("{0}", alglib.ap.format(a2,2)); // EXPECTED: [0,-1,+1] System.Console.ReadLine(); return 0; }
public static int Main(string[] args) { // // Temporaries: // * values of y=x^2-x sampled at three special grids: // * equdistant grid spanning [0,2], x[i] = 2*i/(N-1), i=0..N-1 // * Chebyshev-I grid spanning [-1,+1], x[i] = 1 + Cos(PI*(2*i+1)/(2*n)), i=0..N-1 // * Chebyshev-II grid spanning [-1,+1], x[i] = 1 + Cos(PI*i/(n-1)), i=0..N-1 // * barycentric interpolants for these three grids // * vectors to store coefficients of quadratic representation // double[] y_eqdist = new double[]{0,0,2}; double[] y_cheb1 = new double[]{-0.116025,0.000000,1.616025}; double[] y_cheb2 = new double[]{0,0,2}; alglib.barycentricinterpolant p_eqdist; alglib.barycentricinterpolant p_cheb1; alglib.barycentricinterpolant p_cheb2; double[] a_eqdist; double[] a_cheb1; double[] a_cheb2; // // First, we demonstrate construction of barycentric interpolants on // special grids. We unpack power representation to ensure that // interpolant was built correctly.  // // We calculate interpolant value at x=-2. // In all three cases we should get same f=6 // double t = -2; double v; v = alglib.polynomialcalceqdist(0.0, 2.0, y_eqdist, t); System.Console.WriteLine("{0:F4", v); // EXPECTED: 6.0 v = alglib.polynomialcalccheb1(-1, +1, y_cheb1, t); System.Console.WriteLine("{0:F4", v); // EXPECTED: 6.0 v = alglib.polynomialcalccheb2(-1, +1, y_cheb2, t); System.Console.WriteLine("{0:F4", v); // EXPECTED: 6.0 System.Console.ReadLine(); return 0; }
psif
subpackagepsi
function/************************************************************************* Psi (digamma) function d - psi(x) = -- ln | (x) dx is the logarithmic derivative of the gamma function. For integer x, n-1 - psi(n) = -EUL + > 1/k. - k=1 This formula is used for 0 < n <= 10. If x is negative, it is transformed to a positive argument by the reflection formula psi(1-x) = psi(x) + pi cot(pi x). For general positive x, the argument is made greater than 10 using the recurrence psi(x+1) = psi(x) + 1/x. Then the following asymptotic expansion is applied: inf. B - 2k psi(x) = log(x) - 1/2x - > ------- - 2k k=1 2k x where the B2k are Bernoulli numbers. ACCURACY: Relative error (except absolute when |psi| < 1): arithmetic domain # trials peak rms IEEE 0,30 30000 1.3e-15 1.4e-16 IEEE -30,0 40000 1.5e-15 2.2e-16 Cephes Math Library Release 2.8: June, 2000 Copyright 1984, 1987, 1992, 2000 by Stephen L. Moshier *************************************************************************/public static double alglib.psi(double x)
pspline
subpackagepspline2interpolant
class/************************************************************************* Parametric spline inteprolant: 2-dimensional curve. You should not try to access its members directly - use PSpline2XXXXXXXX() functions instead. *************************************************************************/public class pspline2interpolant { }
pspline3interpolant
class/************************************************************************* Parametric spline inteprolant: 3-dimensional curve. You should not try to access its members directly - use PSpline3XXXXXXXX() functions instead. *************************************************************************/public class pspline3interpolant { }
pspline2arclength
function/************************************************************************* This function calculates arc length, i.e. length of curve between t=a and t=b. INPUT PARAMETERS: P - parametric spline interpolant A,B - parameter values corresponding to arc ends: * B>A will result in positive length returned * B<A will result in negative length returned RESULT: length of arc starting at T=A and ending at T=B. -- ALGLIB PROJECT -- Copyright 30.05.2010 by Bochkanov Sergey *************************************************************************/public static double alglib.pspline2arclength( pspline2interpolant p, double a, double b)
pspline2build
function/************************************************************************* This function builds non-periodic 2-dimensional parametric spline which starts at (X[0],Y[0]) and ends at (X[N-1],Y[N-1]). INPUT PARAMETERS: XY - points, array[0..N-1,0..1]. XY[I,0:1] corresponds to the Ith point. Order of points is important! N - points count, N>=5 for Akima splines, N>=2 for other types of splines. ST - spline type: * 0 Akima spline * 1 parabolically terminated Catmull-Rom spline (Tension=0) * 2 parabolically terminated cubic spline PT - parameterization type: * 0 uniform * 1 chord length * 2 centripetal OUTPUT PARAMETERS: P - parametric spline interpolant NOTES: * this function assumes that there all consequent points are distinct. I.e. (x0,y0)<>(x1,y1), (x1,y1)<>(x2,y2), (x2,y2)<>(x3,y3) and so on. However, non-consequent points may coincide, i.e. we can have (x0,y0)= =(x2,y2). -- ALGLIB PROJECT -- Copyright 28.05.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.pspline2build( double[,] xy, int n, int st, int pt, out pspline2interpolant p)
pspline2buildperiodic
function/************************************************************************* This function builds periodic 2-dimensional parametric spline which starts at (X[0],Y[0]), goes through all points to (X[N-1],Y[N-1]) and then back to (X[0],Y[0]). INPUT PARAMETERS: XY - points, array[0..N-1,0..1]. XY[I,0:1] corresponds to the Ith point. XY[N-1,0:1] must be different from XY[0,0:1]. Order of points is important! N - points count, N>=3 for other types of splines. ST - spline type: * 1 Catmull-Rom spline (Tension=0) with cyclic boundary conditions * 2 cubic spline with cyclic boundary conditions PT - parameterization type: * 0 uniform * 1 chord length * 2 centripetal OUTPUT PARAMETERS: P - parametric spline interpolant NOTES: * this function assumes that there all consequent points are distinct. I.e. (x0,y0)<>(x1,y1), (x1,y1)<>(x2,y2), (x2,y2)<>(x3,y3) and so on. However, non-consequent points may coincide, i.e. we can have (x0,y0)= =(x2,y2). * last point of sequence is NOT equal to the first point. You shouldn't make curve "explicitly periodic" by making them equal. -- ALGLIB PROJECT -- Copyright 28.05.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.pspline2buildperiodic( double[,] xy, int n, int st, int pt, out pspline2interpolant p)
pspline2calc
function/************************************************************************* This function calculates the value of the parametric spline for a given value of parameter T INPUT PARAMETERS: P - parametric spline interpolant T - point: * T in [0,1] corresponds to interval spanned by points * for non-periodic splines T<0 (or T>1) correspond to parts of the curve before the first (after the last) point * for periodic splines T<0 (or T>1) are projected into [0,1] by making T=T-floor(T). OUTPUT PARAMETERS: X - X-position Y - Y-position -- ALGLIB PROJECT -- Copyright 28.05.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.pspline2calc( pspline2interpolant p, double t, out double x, out double y)
pspline2diff
function/************************************************************************* This function calculates derivative, i.e. it returns (dX/dT,dY/dT). INPUT PARAMETERS: P - parametric spline interpolant T - point: * T in [0,1] corresponds to interval spanned by points * for non-periodic splines T<0 (or T>1) correspond to parts of the curve before the first (after the last) point * for periodic splines T<0 (or T>1) are projected into [0,1] by making T=T-floor(T). OUTPUT PARAMETERS: X - X-value DX - X-derivative Y - Y-value DY - Y-derivative -- ALGLIB PROJECT -- Copyright 28.05.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.pspline2diff( pspline2interpolant p, double t, out double x, out double dx, out double y, out double dy)
pspline2diff2
function/************************************************************************* This function calculates first and second derivative with respect to T. INPUT PARAMETERS: P - parametric spline interpolant T - point: * T in [0,1] corresponds to interval spanned by points * for non-periodic splines T<0 (or T>1) correspond to parts of the curve before the first (after the last) point * for periodic splines T<0 (or T>1) are projected into [0,1] by making T=T-floor(T). OUTPUT PARAMETERS: X - X-value DX - derivative D2X - second derivative Y - Y-value DY - derivative D2Y - second derivative -- ALGLIB PROJECT -- Copyright 28.05.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.pspline2diff2( pspline2interpolant p, double t, out double x, out double dx, out double d2x, out double y, out double dy, out double d2y)
pspline2parametervalues
function/************************************************************************* This function returns vector of parameter values correspoding to points. I.e. for P created from (X[0],Y[0])...(X[N-1],Y[N-1]) and U=TValues(P) we have (X[0],Y[0]) = PSpline2Calc(P,U[0]), (X[1],Y[1]) = PSpline2Calc(P,U[1]), (X[2],Y[2]) = PSpline2Calc(P,U[2]), ... INPUT PARAMETERS: P - parametric spline interpolant OUTPUT PARAMETERS: N - array size T - array[0..N-1] NOTES: * for non-periodic splines U[0]=0, U[0]<U[1]<...<U[N-1], U[N-1]=1 * for periodic splines U[0]=0, U[0]<U[1]<...<U[N-1], U[N-1]<1 -- ALGLIB PROJECT -- Copyright 28.05.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.pspline2parametervalues( pspline2interpolant p, out int n, out double[] t)
pspline2tangent
function/************************************************************************* This function calculates tangent vector for a given value of parameter T INPUT PARAMETERS: P - parametric spline interpolant T - point: * T in [0,1] corresponds to interval spanned by points * for non-periodic splines T<0 (or T>1) correspond to parts of the curve before the first (after the last) point * for periodic splines T<0 (or T>1) are projected into [0,1] by making T=T-floor(T). OUTPUT PARAMETERS: X - X-component of tangent vector (normalized) Y - Y-component of tangent vector (normalized) NOTE: X^2+Y^2 is either 1 (for non-zero tangent vector) or 0. -- ALGLIB PROJECT -- Copyright 28.05.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.pspline2tangent( pspline2interpolant p, double t, out double x, out double y)
pspline3arclength
function/************************************************************************* This function calculates arc length, i.e. length of curve between t=a and t=b. INPUT PARAMETERS: P - parametric spline interpolant A,B - parameter values corresponding to arc ends: * B>A will result in positive length returned * B<A will result in negative length returned RESULT: length of arc starting at T=A and ending at T=B. -- ALGLIB PROJECT -- Copyright 30.05.2010 by Bochkanov Sergey *************************************************************************/public static double alglib.pspline3arclength( pspline3interpolant p, double a, double b)
pspline3build
function/************************************************************************* This function builds non-periodic 3-dimensional parametric spline which starts at (X[0],Y[0],Z[0]) and ends at (X[N-1],Y[N-1],Z[N-1]). Same as PSpline2Build() function, but for 3D, so we won't duplicate its description here. -- ALGLIB PROJECT -- Copyright 28.05.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.pspline3build( double[,] xy, int n, int st, int pt, out pspline3interpolant p)
pspline3buildperiodic
function/************************************************************************* This function builds periodic 3-dimensional parametric spline which starts at (X[0],Y[0],Z[0]), goes through all points to (X[N-1],Y[N-1],Z[N-1]) and then back to (X[0],Y[0],Z[0]). Same as PSpline2Build() function, but for 3D, so we won't duplicate its description here. -- ALGLIB PROJECT -- Copyright 28.05.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.pspline3buildperiodic( double[,] xy, int n, int st, int pt, out pspline3interpolant p)
pspline3calc
function/************************************************************************* This function calculates the value of the parametric spline for a given value of parameter T. INPUT PARAMETERS: P - parametric spline interpolant T - point: * T in [0,1] corresponds to interval spanned by points * for non-periodic splines T<0 (or T>1) correspond to parts of the curve before the first (after the last) point * for periodic splines T<0 (or T>1) are projected into [0,1] by making T=T-floor(T). OUTPUT PARAMETERS: X - X-position Y - Y-position Z - Z-position -- ALGLIB PROJECT -- Copyright 28.05.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.pspline3calc( pspline3interpolant p, double t, out double x, out double y, out double z)
pspline3diff
function/************************************************************************* This function calculates derivative, i.e. it returns (dX/dT,dY/dT,dZ/dT). INPUT PARAMETERS: P - parametric spline interpolant T - point: * T in [0,1] corresponds to interval spanned by points * for non-periodic splines T<0 (or T>1) correspond to parts of the curve before the first (after the last) point * for periodic splines T<0 (or T>1) are projected into [0,1] by making T=T-floor(T). OUTPUT PARAMETERS: X - X-value DX - X-derivative Y - Y-value DY - Y-derivative Z - Z-value DZ - Z-derivative -- ALGLIB PROJECT -- Copyright 28.05.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.pspline3diff( pspline3interpolant p, double t, out double x, out double dx, out double y, out double dy, out double z, out double dz)
pspline3diff2
function/************************************************************************* This function calculates first and second derivative with respect to T. INPUT PARAMETERS: P - parametric spline interpolant T - point: * T in [0,1] corresponds to interval spanned by points * for non-periodic splines T<0 (or T>1) correspond to parts of the curve before the first (after the last) point * for periodic splines T<0 (or T>1) are projected into [0,1] by making T=T-floor(T). OUTPUT PARAMETERS: X - X-value DX - derivative D2X - second derivative Y - Y-value DY - derivative D2Y - second derivative Z - Z-value DZ - derivative D2Z - second derivative -- ALGLIB PROJECT -- Copyright 28.05.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.pspline3diff2( pspline3interpolant p, double t, out double x, out double dx, out double d2x, out double y, out double dy, out double d2y, out double z, out double dz, out double d2z)
pspline3parametervalues
function/************************************************************************* This function returns vector of parameter values correspoding to points. Same as PSpline2ParameterValues(), but for 3D. -- ALGLIB PROJECT -- Copyright 28.05.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.pspline3parametervalues( pspline3interpolant p, out int n, out double[] t)
pspline3tangent
function/************************************************************************* This function calculates tangent vector for a given value of parameter T INPUT PARAMETERS: P - parametric spline interpolant T - point: * T in [0,1] corresponds to interval spanned by points * for non-periodic splines T<0 (or T>1) correspond to parts of the curve before the first (after the last) point * for periodic splines T<0 (or T>1) are projected into [0,1] by making T=T-floor(T). OUTPUT PARAMETERS: X - X-component of tangent vector (normalized) Y - Y-component of tangent vector (normalized) Z - Z-component of tangent vector (normalized) NOTE: X^2+Y^2+Z^2 is either 1 (for non-zero tangent vector) or 0. -- ALGLIB PROJECT -- Copyright 28.05.2010 by Bochkanov Sergey *************************************************************************/public static void alglib.pspline3tangent( pspline3interpolant p, double t, out double x, out double y, out double z)
ratint
subpackagebarycentricinterpolant
class/************************************************************************* Barycentric interpolant. *************************************************************************/public class barycentricinterpolant { }
barycentricbuildfloaterhormann
function/************************************************************************* Rational interpolant without poles The subroutine constructs the rational interpolating function without real poles (see 'Barycentric rational interpolation with no poles and high rates of approximation', Michael S. Floater. and Kai Hormann, for more information on this subject). Input parameters: X - interpolation nodes, array[0..N-1]. Y - function values, array[0..N-1]. N - number of nodes, N>0. D - order of the interpolation scheme, 0 <= D <= N-1. D<0 will cause an error. D>=N it will be replaced with D=N-1. if you don't know what D to choose, use small value about 3-5. Output parameters: B - barycentric interpolant. Note: this algorithm always succeeds and calculates the weights with close to machine precision. -- ALGLIB PROJECT -- Copyright 17.06.2007 by Bochkanov Sergey *************************************************************************/public static void alglib.barycentricbuildfloaterhormann( double[] x, double[] y, int n, int d, out barycentricinterpolant b)
barycentricbuildxyw
function/************************************************************************* Rational interpolant from X/Y/W arrays F(t) = SUM(i=0,n-1,w[i]*f[i]/(t-x[i])) / SUM(i=0,n-1,w[i]/(t-x[i])) INPUT PARAMETERS: X - interpolation nodes, array[0..N-1] F - function values, array[0..N-1] W - barycentric weights, array[0..N-1] N - nodes count, N>0 OUTPUT PARAMETERS: B - barycentric interpolant built from (X, Y, W) -- ALGLIB -- Copyright 17.08.2009 by Bochkanov Sergey *************************************************************************/public static void alglib.barycentricbuildxyw( double[] x, double[] y, double[] w, int n, out barycentricinterpolant b)
barycentriccalc
function/************************************************************************* Rational interpolation using barycentric formula F(t) = SUM(i=0,n-1,w[i]*f[i]/(t-x[i])) / SUM(i=0,n-1,w[i]/(t-x[i])) Input parameters: B - barycentric interpolant built with one of model building subroutines. T - interpolation point Result: barycentric interpolant F(t) -- ALGLIB -- Copyright 17.08.2009 by Bochkanov Sergey *************************************************************************/public static double alglib.barycentriccalc( barycentricinterpolant b, double t)
barycentricdiff1
function/************************************************************************* Differentiation of barycentric interpolant: first derivative. Algorithm used in this subroutine is very robust and should not fail until provided with values too close to MaxRealNumber (usually MaxRealNumber/N or greater will overflow). INPUT PARAMETERS: B - barycentric interpolant built with one of model building subroutines. T - interpolation point OUTPUT PARAMETERS: F - barycentric interpolant at T DF - first derivative NOTE -- ALGLIB -- Copyright 17.08.2009 by Bochkanov Sergey *************************************************************************/public static void alglib.barycentricdiff1( barycentricinterpolant b, double t, out double f, out double df)
barycentricdiff2
function/************************************************************************* Differentiation of barycentric interpolant: first/second derivatives. INPUT PARAMETERS: B - barycentric interpolant built with one of model building subroutines. T - interpolation point OUTPUT PARAMETERS: F - barycentric interpolant at T DF - first derivative D2F - second derivative NOTE: this algorithm may fail due to overflow/underflor if used on data whose values are close to MaxRealNumber or MinRealNumber. Use more robust BarycentricDiff1() subroutine in such cases. -- ALGLIB -- Copyright 17.08.2009 by Bochkanov Sergey *************************************************************************/public static void alglib.barycentricdiff2( barycentricinterpolant b, double t, out double f, out double df, out double d2f)
barycentriclintransx
function/************************************************************************* This subroutine performs linear transformation of the argument. INPUT PARAMETERS: B - rational interpolant in barycentric form CA, CB - transformation coefficients: x = CA*t + CB OUTPUT PARAMETERS: B - transformed interpolant with X replaced by T -- ALGLIB PROJECT -- Copyright 19.08.2009 by Bochkanov Sergey *************************************************************************/public static void alglib.barycentriclintransx( barycentricinterpolant b, double ca, double cb)
barycentriclintransy
function/************************************************************************* This subroutine performs linear transformation of the barycentric interpolant. INPUT PARAMETERS: B - rational interpolant in barycentric form CA, CB - transformation coefficients: B2(x) = CA*B(x) + CB OUTPUT PARAMETERS: B - transformed interpolant -- ALGLIB PROJECT -- Copyright 19.08.2009 by Bochkanov Sergey *************************************************************************/public static void alglib.barycentriclintransy( barycentricinterpolant b, double ca, double cb)
barycentricunpack
function/************************************************************************* Extracts X/Y/W arrays from rational interpolant INPUT PARAMETERS: B - barycentric interpolant OUTPUT PARAMETERS: N - nodes count, N>0 X - interpolation nodes, array[0..N-1] F - function values, array[0..N-1] W - barycentric weights, array[0..N-1] -- ALGLIB -- Copyright 17.08.2009 by Bochkanov Sergey *************************************************************************/public static void alglib.barycentricunpack( barycentricinterpolant b, out int n, out double[] x, out double[] y, out double[] w)
rbf
subpackage| rbf_d_ml_ls | Least squares problem solved with RBF-ML algorithm | |
| rbf_d_ml_simple | Simple model built with RBF-ML algorithm | |
| rbf_d_polterm | RBF models - working with polynomial term | |
| rbf_d_qnn | Simple model built with RBF-QNN algorithm | |
| rbf_d_serialize | Serialization/unserialization | |
| rbf_d_vector | Working with vector functions |
rbfmodel
class/************************************************************************* RBF model. Never try to directly work with fields of this object - always use ALGLIB functions to use this object. *************************************************************************/public class rbfmodel { }
rbfreport
class/************************************************************************* RBF solution report: * TerminationType - termination type, positive values - success, non-positive - failure. *************************************************************************/public class rbfreport { public int arows; public int acols; public int annz; public int iterationscount; public int nmv; public int terminationtype; }
rbfbuildmodel
function/************************************************************************* This function builds RBF model and returns report (contains some information which can be used for evaluation of the algorithm properties). Call to this function modifies RBF model by calculating its centers/radii/ weights and saving them into RBFModel structure. Initially RBFModel contain zero coefficients, but after call to this function we will have coefficients which were calculated in order to fit our dataset. After you called this function you can call RBFCalc(), RBFGridCalc() and other model calculation functions. INPUT PARAMETERS: S - RBF model, initialized by RBFCreate() call Rep - report: * Rep.TerminationType: * -5 - non-distinct basis function centers were detected, interpolation aborted * -4 - nonconvergence of the internal SVD solver * 1 - successful termination Fields are used for debugging purposes: * Rep.IterationsCount - iterations count of the LSQR solver * Rep.NMV - number of matrix-vector products * Rep.ARows - rows count for the system matrix * Rep.ACols - columns count for the system matrix * Rep.ANNZ - number of significantly non-zero elements (elements above some algorithm-determined threshold) NOTE: failure to build model will leave current state of the structure unchanged. -- ALGLIB -- Copyright 13.12.2011 by Bochkanov Sergey *************************************************************************/public static void alglib.rbfbuildmodel(rbfmodel s, out rbfreport rep)
rbfcalc
function/************************************************************************* This function calculates values of the RBF model at the given point. This is general function which can be used for arbitrary NX (dimension of the space of arguments) and NY (dimension of the function itself). However when you have NY=1 you may find more convenient to use RBFCalc2() or RBFCalc3(). This function returns 0.0 when model is not initialized. INPUT PARAMETERS: S - RBF model X - coordinates, array[NX]. X may have more than NX elements, in this case only leading NX will be used. OUTPUT PARAMETERS: Y - function value, array[NY]. Y is out-parameter and reallocated after call to this function. In case you want to reuse previously allocated Y, you may use RBFCalcBuf(), which reallocates Y only when it is too small. -- ALGLIB -- Copyright 13.12.2011 by Bochkanov Sergey *************************************************************************/public static void alglib.rbfcalc(rbfmodel s, double[] x, out double[] y)
Examples: [1]
rbfcalc2
function/************************************************************************* This function calculates values of the RBF model in the given point. This function should be used when we have NY=1 (scalar function) and NX=2 (2-dimensional space). If you have 3-dimensional space, use RBFCalc3(). If you have general situation (NX-dimensional space, NY-dimensional function) you should use general, less efficient implementation RBFCalc(). If you want to calculate function values many times, consider using RBFGridCalc2(), which is far more efficient than many subsequent calls to RBFCalc2(). This function returns 0.0 when: * model is not initialized * NX<>2 *NY<>1 INPUT PARAMETERS: S - RBF model X0 - first coordinate, finite number X1 - second coordinate, finite number RESULT: value of the model or 0.0 (as defined above) -- ALGLIB -- Copyright 13.12.2011 by Bochkanov Sergey *************************************************************************/public static double alglib.rbfcalc2(rbfmodel s, double x0, double x1)
rbfcalc3
function/************************************************************************* This function calculates values of the RBF model in the given point. This function should be used when we have NY=1 (scalar function) and NX=3 (3-dimensional space). If you have 2-dimensional space, use RBFCalc2(). If you have general situation (NX-dimensional space, NY-dimensional function) you should use general, less efficient implementation RBFCalc(). This function returns 0.0 when: * model is not initialized * NX<>3 *NY<>1 INPUT PARAMETERS: S - RBF model X0 - first coordinate, finite number X1 - second coordinate, finite number X2 - third coordinate, finite number RESULT: value of the model or 0.0 (as defined above) -- ALGLIB -- Copyright 13.12.2011 by Bochkanov Sergey *************************************************************************/public static double alglib.rbfcalc3( rbfmodel s, double x0, double x1, double x2)
rbfcalcbuf
function/************************************************************************* This function calculates values of the RBF model at the given point. Same as RBFCalc(), but does not reallocate Y when in is large enough to store function values. INPUT PARAMETERS: S - RBF model X - coordinates, array[NX]. X may have more than NX elements, in this case only leading NX will be used. Y - possibly preallocated array OUTPUT PARAMETERS: Y - function value, array[NY]. Y is not reallocated when it is larger than NY. -- ALGLIB -- Copyright 13.12.2011 by Bochkanov Sergey *************************************************************************/public static void alglib.rbfcalcbuf( rbfmodel s, double[] x, ref double[] y)
rbfcreate
function/************************************************************************* This function creates RBF model for a scalar (NY=1) or vector (NY>1) function in a NX-dimensional space (NX=2 or NX=3). Newly created model is empty. It can be used for interpolation right after creation, but it just returns zeros. You have to add points to the model, tune interpolation settings, and then call model construction function RBFBuildModel() which will update model according to your specification. USAGE: 1. User creates model with RBFCreate() 2. User adds dataset with RBFSetPoints() (points do NOT have to be on a regular grid) 3. (OPTIONAL) User chooses polynomial term by calling: * RBFLinTerm() to set linear term * RBFConstTerm() to set constant term * RBFZeroTerm() to set zero term By default, linear term is used. 4. User chooses specific RBF algorithm to use: either QNN (RBFSetAlgoQNN) or ML (RBFSetAlgoMultiLayer). 5. User calls RBFBuildModel() function which rebuilds model according to the specification 6. User may call RBFCalc() to calculate model value at the specified point, RBFGridCalc() to calculate model values at the points of the regular grid. User may extract model coefficients with RBFUnpack() call. INPUT PARAMETERS: NX - dimension of the space, NX=2 or NX=3 NY - function dimension, NY>=1 OUTPUT PARAMETERS: S - RBF model (initially equals to zero) NOTE 1: memory requirements. RBF models require amount of memory which is proportional to the number of data points. Memory is allocated during model construction, but most of this memory is freed after model coefficients are calculated. Some approximate estimates for N centers with default settings are given below: * about 250*N*(sizeof(double)+2*sizeof(int)) bytes of memory is needed during model construction stage. * about 15*N*sizeof(double) bytes is needed after model is built. For example, for N=100000 we may need 0.6 GB of memory to build model, but just about 0.012 GB to store it. -- ALGLIB -- Copyright 13.12.2011 by Bochkanov Sergey *************************************************************************/public static void alglib.rbfcreate(int nx, int ny, out rbfmodel s)
Examples: [1] [2] [3] [4] [5] [6]
rbfgridcalc2
function/************************************************************************* This function calculates values of the RBF model at the regular grid. Grid have N0*N1 points, with Point[I,J] = (X0[I], X1[J]) This function returns 0.0 when: * model is not initialized * NX<>2 *NY<>1 INPUT PARAMETERS: S - RBF model X0 - array of grid nodes, first coordinates, array[N0] N0 - grid size (number of nodes) in the first dimension X1 - array of grid nodes, second coordinates, array[N1] N1 - grid size (number of nodes) in the second dimension OUTPUT PARAMETERS: Y - function values, array[N0,N1]. Y is out-variable and is reallocated by this function. -- ALGLIB -- Copyright 13.12.2011 by Bochkanov Sergey *************************************************************************/public static void alglib.rbfgridcalc2( rbfmodel s, double[] x0, int n0, double[] x1, int n1, out double[,] y)
rbfserialize
function/************************************************************************* This function serializes data structure to string. Important properties of s_out: * it contains alphanumeric characters, dots, underscores, minus signs * these symbols are grouped into words, which are separated by spaces and Windows-style (CR+LF) newlines * although serializer uses spaces and CR+LF as separators, you can replace any separator character by arbitrary combination of spaces, tabs, Windows or Unix newlines. It allows flexible reformatting of the string in case you want to include it into text or XML file. But you should not insert separators into the middle of the "words" nor you should change case of letters. * s_out can be freely moved between 32-bit and 64-bit systems, little and big endian machines, and so on. You can serialize structure on 32-bit machine and unserialize it on 64-bit one (or vice versa), or serialize it on SPARC and unserialize on x86. You can also serialize it in C++ version of ALGLIB and unserialize in C# one, and vice versa. *************************************************************************/public static void rbfserialize(rbfmodel obj, out string s_out)
rbfsetalgomultilayer
function/************************************************************************* This function sets RBF interpolation algorithm. ALGLIB supports several RBF algorithms with different properties. This algorithm is called RBF-ML. It builds multilayer RBF model, i.e. model with subsequently decreasing radii, which allows us to combine smoothness (due to large radii of the first layers) with exactness (due to small radii of the last layers) and fast convergence. Internally RBF-ML uses many different means of acceleration, from sparse matrices to KD-trees, which results in algorithm whose working time is roughly proportional to N*log(N)*Density*RBase^2*NLayers, where N is a number of points, Density is an average density if points per unit of the interpolation space, RBase is an initial radius, NLayers is a number of layers. RBF-ML is good for following kinds of interpolation problems: 1. "exact" problems (perfect fit) with well separated points 2. least squares problems with arbitrary distribution of points (algorithm gives perfect fit where it is possible, and resorts to least squares fit in the hard areas). 3. noisy problems where we want to apply some controlled amount of smoothing. INPUT PARAMETERS: S - RBF model, initialized by RBFCreate() call RBase - RBase parameter, RBase>0 NLayers - NLayers parameter, NLayers>0, recommended value to start with - about 5. LambdaV - regularization value, can be useful when solving problem in the least squares sense. Optimal lambda is problem- dependent and require trial and error. In our experience, good lambda can be as large as 0.1, and you can use 0.001 as initial guess. Default value - 0.01, which is used when LambdaV is not given. You can specify zero value, but it is not recommended to do so. TUNING ALGORITHM In order to use this algorithm you have to choose three parameters: * initial radius RBase * number of layers in the model NLayers * regularization coefficient LambdaV Initial radius is easy to choose - you can pick any number several times larger than the average distance between points. Algorithm won't break down if you choose radius which is too large (model construction time will increase, but model will be built correctly). Choose such number of layers that RLast=RBase/2^(NLayers-1) (radius used by the last layer) will be smaller than the typical distance between points. In case model error is too large, you can increase number of layers. Having more layers will make model construction and evaluation proportionally slower, but it will allow you to have model which precisely fits your data. From the other side, if you want to suppress noise, you can DECREASE number of layers to make your model less flexible. Regularization coefficient LambdaV controls smoothness of the individual models built for each layer. We recommend you to use default value in case you don't want to tune this parameter, because having non-zero LambdaV accelerates and stabilizes internal iterative algorithm. In case you want to suppress noise you can use LambdaV as additional parameter (larger value = more smoothness) to tune. TYPICAL ERRORS 1. Using initial radius which is too large. Memory requirements of the RBF-ML are roughly proportional to N*Density*RBase^2 (where Density is an average density of points per unit of the interpolation space). In the extreme case of the very large RBase we will need O(N^2) units of memory - and many layers in order to decrease radius to some reasonably small value. 2. Using too small number of layers - RBF models with large radius are not flexible enough to reproduce small variations in the target function. You need many layers with different radii, from large to small, in order to have good model. 3. Using initial radius which is too small. You will get model with "holes" in the areas which are too far away from interpolation centers. However, algorithm will work correctly (and quickly) in this case. 4.  This model will perfectly reproduce your function, but maybe you will be able to achieve similar results with less layers (and less memory). -- ALGLIB -- Copyright 02.03.2012 by Bochkanov Sergey *************************************************************************/public static void alglib.rbfsetalgomultilayer( rbfmodel s, double rbase, int nlayers) public static void alglib.rbfsetalgomultilayer( rbfmodel s, double rbase, int nlayers, double lambdav)
rbfsetalgoqnn
function/************************************************************************* This function sets RBF interpolation algorithm. ALGLIB supports several RBF algorithms with different properties. This algorithm is called RBF-QNN and it is good for point sets with following properties: a) all points are distinct b) all points are well separated. c) points distribution is approximately uniform. There is no "contour lines", clusters of points, or other small-scale structures. Algorithm description: 1) interpolation centers are allocated to data points 2) interpolation radii are calculated as distances to the nearest centers times Q coefficient (where Q is a value from [0.75,1.50]). 3) after performing (2) radii are transformed in order to avoid situation when single outlier has very large radius and influences many points across all dataset. Transformation has following form: new_r[i] = min(r[i],Z*median(r[])) where r[i] is I-th radius, median() is a median radius across entire dataset, Z is user-specified value which controls amount of deviation from median radius. When (a) is violated, we will be unable to build RBF model. When (b) or (c) are violated, model will be built, but interpolation quality will be low. See http://www.alglib.net/interpolation/ for more information on this subject. This algorithm is used by default. Additional Q parameter controls smoothness properties of the RBF basis: * Q<0.75 will give perfectly conditioned basis, but terrible smoothness properties (RBF interpolant will have sharp peaks around function values) * Q around 1.0 gives good balance between smoothness and condition number * Q>1.5 will lead to badly conditioned systems and slow convergence of the underlying linear solver (although smoothness will be very good) * Q>2.0 will effectively make optimizer useless because it won't converge within reasonable amount of iterations. It is possible to set such large Q, but it is advised not to do so. INPUT PARAMETERS: S - RBF model, initialized by RBFCreate() call Q - Q parameter, Q>0, recommended value - 1.0 Z - Z parameter, Z>0, recommended value - 5.0 NOTE: this function has some serialization-related subtleties. We recommend you to study serialization examples from ALGLIB Reference Manual if you want to perform serialization of your models. -- ALGLIB -- Copyright 13.12.2011 by Bochkanov Sergey *************************************************************************/public static void alglib.rbfsetalgoqnn(rbfmodel s) public static void alglib.rbfsetalgoqnn(rbfmodel s, double q, double z)
rbfsetconstterm
function/************************************************************************* This function sets constant term (model is a sum of radial basis functions plus constant). This function won't have effect until next call to RBFBuildModel(). INPUT PARAMETERS: S - RBF model, initialized by RBFCreate() call NOTE: this function has some serialization-related subtleties. We recommend you to study serialization examples from ALGLIB Reference Manual if you want to perform serialization of your models. -- ALGLIB -- Copyright 13.12.2011 by Bochkanov Sergey *************************************************************************/public static void alglib.rbfsetconstterm(rbfmodel s)
Examples: [1]
rbfsetlinterm
function/************************************************************************* This function sets linear term (model is a sum of radial basis functions plus linear polynomial). This function won't have effect until next call to RBFBuildModel(). INPUT PARAMETERS: S - RBF model, initialized by RBFCreate() call NOTE: this function has some serialization-related subtleties. We recommend you to study serialization examples from ALGLIB Reference Manual if you want to perform serialization of your models. -- ALGLIB -- Copyright 13.12.2011 by Bochkanov Sergey *************************************************************************/public static void alglib.rbfsetlinterm(rbfmodel s)
Examples: [1]
rbfsetpoints
function/************************************************************************* This function adds dataset. This function overrides results of the previous calls, i.e. multiple calls of this function will result in only the last set being added. INPUT PARAMETERS: S - RBF model, initialized by RBFCreate() call. XY - points, array[N,NX+NY]. One row corresponds to one point in the dataset. First NX elements are coordinates, next NY elements are function values. Array may be larger than specific, in this case only leading [N,NX+NY] elements will be used. N - number of points in the dataset After you've added dataset and (optionally) tuned algorithm settings you should call RBFBuildModel() in order to build a model for you. NOTE: this function has some serialization-related subtleties. We recommend you to study serialization examples from ALGLIB Reference Manual if you want to perform serialization of your models. -- ALGLIB -- Copyright 13.12.2011 by Bochkanov Sergey *************************************************************************/public static void alglib.rbfsetpoints(rbfmodel s, double[,] xy) public static void alglib.rbfsetpoints(rbfmodel s, double[,] xy, int n)
rbfsetzeroterm
function/************************************************************************* This function sets zero term (model is a sum of radial basis functions without polynomial term). This function won't have effect until next call to RBFBuildModel(). INPUT PARAMETERS: S - RBF model, initialized by RBFCreate() call NOTE: this function has some serialization-related subtleties. We recommend you to study serialization examples from ALGLIB Reference Manual if you want to perform serialization of your models. -- ALGLIB -- Copyright 13.12.2011 by Bochkanov Sergey *************************************************************************/public static void alglib.rbfsetzeroterm(rbfmodel s)
Examples: [1]
rbfunpack
function/************************************************************************* This function "unpacks" RBF model by extracting its coefficients. INPUT PARAMETERS: S - RBF model OUTPUT PARAMETERS: NX - dimensionality of argument NY - dimensionality of the target function XWR - model information, array[NC,NX+NY+1]. One row of the array corresponds to one basis function: * first NX columns - coordinates of the center * next NY columns - weights, one per dimension of the function being modelled * last column - radius, same for all dimensions of the function being modelled NC - number of the centers V - polynomial term , array[NY,NX+1]. One row per one dimension of the function being modelled. First NX elements are linear coefficients, V[NX] is equal to the constant part. -- ALGLIB -- Copyright 13.12.2011 by Bochkanov Sergey *************************************************************************/public static void alglib.rbfunpack( rbfmodel s, out int nx, out int ny, out double[,] xwr, out int nc, out double[,] v)
Examples: [1]
rbfunserialize
function/************************************************************************* This function unserializes data structure from string. *************************************************************************/public static void rbfunserialize(string s_in, out rbfmodel obj)
public static int Main(string[] args) { // // This example shows how to solve least squares problems with RBF-ML algorithm. // Below we assume that you already know basic concepts shown in the RBF_D_QNN and // RBF_D_ML_SIMPLE examples. // alglib.rbfmodel model; alglib.rbfreport rep; double v; // // We have 2-dimensional space and very simple fitting problem - all points // except for two are well separated and located at straight line. Two // "exceptional" points are very close, with distance between them as small // as 0.01. RBF-QNN algorithm will have many difficulties with such distribution // of points: // X Y // -2 1 // -1 0 // -0.005 1 // +0.005 2 // +1 -1 // +2 1 // How will RBF-ML handle such problem? // alglib.rbfcreate(2, 1, out model); double[,] xy0 = new double[,]{{-2,0,1},{-1,0,0},{-0.005,0,1},{+0.005,0,2},{+1,0,-1},{+2,0,1}}; alglib.rbfsetpoints(model, xy0); // First, we try to use R=5.0 with single layer (NLayers=1) and moderate amount // of regularization. Well, we already expected that results will be bad: // Model(x=-2,y=0)=0.8407 (instead of 1.0) // Model(x=0.005,y=0)=0.6584 (instead of 2.0) // We need more layers to show better results. alglib.rbfsetalgomultilayer(model, 5.0, 1, 1.0e-3); alglib.rbfbuildmodel(model, out rep); v = alglib.rbfcalc2(model, -2.0, 0.0); System.Console.WriteLine("{0:F2", v); // EXPECTED: 0.8407981659 v = alglib.rbfcalc2(model, 0.005, 0.0); System.Console.WriteLine("{0:F2", v); // EXPECTED: 0.6584267649 // With 4 layers we got better result at x=-2 (point which is well separated // from its neighbors). Model is now many times closer to the original data // Model(x=-2,y=0)=0.9992 (instead of 1.0) // Model(x=0.005,y=0)=1.5534 (instead of 2.0) // We may see that at x=0.005 result is a bit closer to 2.0, but does not // reproduce function value precisely because of close neighbor located at // at x=-0.005. Let's add two layers... alglib.rbfsetalgomultilayer(model, 5.0, 4, 1.0e-3); alglib.rbfbuildmodel(model, out rep); v = alglib.rbfcalc2(model, -2.0, 0.0); System.Console.WriteLine("{0:F2", v); // EXPECTED: 0.9992673278 v = alglib.rbfcalc2(model, 0.005, 0.0); System.Console.WriteLine("{0:F2", v); // EXPECTED: 1.5534666012 // With 6 layers we got almost perfect fit: // Model(x=-2,y=0)=1.000 (perfect fit) // Model(x=0.005,y=0)=1.996 (instead of 2.0) // Of course, we can reduce error at x=0.005 down to zero by adding more // layers. But do we really need it? alglib.rbfsetalgomultilayer(model, 5.0, 6, 1.0e-3); alglib.rbfbuildmodel(model, out rep); v = alglib.rbfcalc2(model, -2.0, 0.0); System.Console.WriteLine("{0:F2", v); // EXPECTED: 1.0000000000 v = alglib.rbfcalc2(model, 0.005, 0.0); System.Console.WriteLine("{0:F2", v); // EXPECTED: 1.9965775952 // Do we really need zero error? We have f(+0.005)=2 and f(-0.005)=1. // Two points are very close, and in real life situations it often means // that difference in function values can be explained by noise in the // data. Thus, true value of the underlying function should be close to // 1.5 (halfway between 1.0 and 2.0). // // How can we get such result with RBF-ML? Well, we can: // a) reduce number of layers (make model less flexible) // b) increase regularization coefficient (another way of reducing flexibility) // // Having NLayers=5 and LambdaV=0.1 gives us good least squares fit to the data: // Model(x=-2,y=0)=1.000 // Model(x=-0.005,y=0)=1.504 // Model(x=+0.005,y=0)=1.496 alglib.rbfsetalgomultilayer(model, 5.0, 5, 1.0e-1); alglib.rbfbuildmodel(model, out rep); v = alglib.rbfcalc2(model, -2.0, 0.0); System.Console.WriteLine("{0:F2", v); // EXPECTED: 1.0000001620 v = alglib.rbfcalc2(model, -0.005, 0.0); System.Console.WriteLine("{0:F2", v); // EXPECTED: 1.5042954378 v = alglib.rbfcalc2(model, 0.005, 0.0); System.Console.WriteLine("{0:F2", v); // EXPECTED: 1.4957042013 System.Console.ReadLine(); return 0; }