|
||||||||
| PREV NEXT | FRAMES NO FRAMES | |||||||
See:
Description
| Core Packages | |
| ru.sscc.util | Collection of miscellaneous utility classes. |
| ru.sscc.util.data | Collection of data access classes implementing such concepts as a data container, data vector, data pointer, and array of vectors. |
| Linear Algebra | |
| ru.sscc.matrix | Collection of classes describing algebraic matrices. |
| ru.sscc.matrix.solve | Collection of Solvers intended for solving of Systems of Linear Algebraic Equations (SLAE). |
| Splines | |
| ru.sscc.spline | Basic spline approximation classes. |
| ru.sscc.spline.analytic | Multivariate Duchon's splines. |
| ru.sscc.spline.creator | The collection of auxiliary classes used in spline creation. |
| ru.sscc.spline.mesh | The collection of 1D meshes and auxiliary classes. |
| ru.sscc.spline.polynomial | One-dimensional polynomial splines. |
| ru.sscc.spline.reduction | The collection of multi-D space transformations (reductions) and multi-D meshes. |
Version 1.25.
See also: JSpline+: Change List.
The samples and benchmarks are described in Samples for JSpline+
The first and main aim of the development of this Library was an implementation of various splines to be constructed by function measurements on 1D or multi-dimensional mesh of scattered points. This job was originated on the spline software library LIDA-3 developed in 1985 at the Novosibirsk Computing Center under the leadership of Prof. Vladimir A. Vasilenko. The current job began at the middle of 1999. While developing, the initial scope of the Library was enlarged. So, the "plus" sign in the library name reflects the fact that the library contains other independent sections within. It now consists of 3 main sections:
The scope of the Library will be continuously enlarged to cover other fields of the Numerical Mathematics.
The ru.sscc.util.data package provides unified access to user's data vectors. The enumeration of vector's entries always starts from zero. All package classes support four primitive operations upon data entries:
Four conceptual classes are introduced here:
RealContainerrelativeAccuracy()
method returns the relative accuracy of data type used in the
container, e.g. the smallest value e such that 1+e
isn't equal to 1 for the data type used. This constant is widely useful in
numeric calculations, but it is absent in Java.
RealVector
RealPointerSimplePointer,
and to point to entries of the RealVectors instance, the
RealPointers.
RealVectorssquaredDistance(index,pointer)
method is used.
All above mentioned classes are tied together: the RealContainer class allows to create a vector or pointer upon it; the RealVector class allows to create a subvector or pointer to a subvector in it; the RealVectors class provides getting of i-th vector from the set of vectors, creation of a vector consisting of j-th component of all vectors, creation of a pointer to i-th vector. All (sub)vector creation operations don't create a copy of entries. They create an instance that refers to the same underlying data container as the instance with the help of which the operation was performed.
Two special subclasses,
ChainVector and
ChainPointer, provide the access to the chain of
vectors considered as one composite vector.
A container can move or shift the starting index of a vector created by it. This operation is useful in matrix calculations when a matrix is considered as a set of rows or columns, and the same operation must be produced on every element of the set.
To simplify conversions, a number of static by methods is introduced.
For example, to create a RealContainer from a double[] or
float[] array, you can write down the following:
double[] array = new double[count];
//... Fill the array in
RealContainer container = RealContainer.by(array);
Finally, we notice the RealMath class. It
contains many useful static methods that perform some BLAS operations upon
double[], float[], and RealPointer instances. The
calculation of factorial of n, integer power of x, and a polynomial value
by the Gorner's algorithm is also implemented here.
The ru.sscc.matrix package and its subpackages present the Numerical Linear Algebra Section of the Library. The root package describes matrices. Only real matrices are implemented now.
The Matrix interface describes the most common properties
of any matrix: a matrix is a rectangular table consisting of a constant number
of rows and columns and a matrix is serializable object.
The RealAlgebraicMatrix interface extends
the Matrix by adding two algebraic methods: multiplication of the matrix
by a vector and multiplication of the transposed matrix by a vector. In other
words, the algebraic matrix is the operator that linearly maps a source vector
to the target vector and it has the transposed operator acting analogously. In
addition, the method relativeAccuracy() must return the relative
accuracy of data type used for matrix entries.
When a matrix is used for solving of a linear algebraic system, the decomposition coefficients are stored in it. In this case, the matrix loses the algebraicity and the multiplication by a vector has no sense more. We say that the matrix becomes non-algebraic and the try to multiply it by a vector throws the IllegalStateException.
The TransposedMatrix class is the decorator
class. It is used if you need to "transpose" an algebraic matrix. This is the
logical transposing. It logically interchanges the columns and rows numbers and
treats the matrix-by-vector multiplication as the transposed multiplication for
the underlying matrix. Conversely, the multiplication of the transposed matrix
by a vector is treated as the direct multiplication for the underlying matrix.
The RealMatrix class is the basic class for almost
all real matrices. It is the abstract class. It partially implements the
RealAlgebraicMatrix interface and adds four primitive operations for the
access to matrix values,
get(i,j),
set(i,j,value),
add(i,j,value), and
mul(i,j,value).
The indices of matrix entries start from zero. So, if the matrix has m rows and n columns then i must belong to the set {0,1,...,m-1} and j must belong to the set {0,1,...,n-1}.
All real matrices store their coefficients in the body, an instance of
the
RealContainer type, e.g. the 1D array is used as
the storage for matrix entries. This solution was selected by the efficiency
reason. You can get the matrix body by the
getContainer()
method. The access to the matrix dimensions may be done in two way:
rowsNumber() and
columnsNumber() methods inherited
from the Matrix interface and
nRows and
nColumns.
The RealMatrix class supports cloning. The clone is the matrix having the separate body copied from the body of the original matrix. The clone is the algebraic matrix, independently to the state of the original matrix. If the original matrix is a submatrix of a dense matrix, the clone will contain the submatrix entries only.
Three types of matrices are implemented in the current version of the Library:
DenseMatrix
SymBandedMatrix
RectBandedMatrixThe switching between algebraic and non-algebraic state of matrices is implemented by the lock and reuse methods: the lock method changes the matrix state to be non-algebraic and the reuse method allows the reusing of the matrix as the algebraic matrix.
The package supports matrix-by-matrix multiplication operations as static
methods of the RealMatrix class.
The ru.sscc.matrix.solve package is devoted to the solving of a System of Linear Algebraic Equations (SLAE)
|
Every real solver implements the
RealSolver
interface, containing the
solve(b,x)
method that receives the right hand side vector b and calculates the solution
vector x. Two auxiliary methods are also described in the interface:
sourceSize() and
targetSize().
The RealSolver interfaces supports serialization.
We distinct two main types of real solvers: direct solvers and iterative
solvers. Only direct solvers are implemented yet. The basic abstract class for
all direct solvers is the
RealDirectSolver.
It implements the basic interface for all direct solvers and adds the abstract
method
factorize()
that provides the matrix decomposition by factors (e.g. factorization). The
matrix to be factorized is attached to the solver. A direct solver can have two
states: before the factorization and after it. If the matrix is nonfactorized
yet, the solve method is inaccessible. When the matrix is factorized, it
becomes non-algebraic, and all algebraic operations with it are blocked. The
reuse method implemented in all subclasses changes the state of the
solver and the attached matrix to the initial state: the factorization tag is
cleared and the matrix becomes algebraic again. So, after that, you can fill in
the attached matrix again.
The basic principal for all direct solvers is: "Once factorize, many times solve". So, if you need to solve many SLAEs that differ in the right hand sides only, you should factorize the matrix only once, and then can find the solution many times for different right hand sides.
To improve the solution accuracy, use the iterative refinement algorithm implemented in the RealDirectSolver. It is most efficient if the attached matrix is constructed on the base of the float[] type vector. The construction of the inverse matrix is also implemented in this class.
The RealDirectSolver class supports cloning. A clone is always clear solver without attached matrix.
The hierarchy of the classes based on the RealDirectSolver is the following:
CholeskyBandSolver
that solves the SLAE with a positive definite symmetric banded matrix by
Cholesky Method also known as Square Root Method.
RealDenseSolver
is the basic class for solvers oriented to the systems with dense matrices. It
contains many subclasses. The direct non-abstract subclasses intended for the
solving of SLAE with a square dense matrix are: (1) the
CholeskySolver
solves the system with a positive definite symmetric matrix; (2) the
GaussSolver
solves the system with a positive or negative definite non-symmetric matrix
(in this case the Gauss eliminations may be produced directly and pivoting
isn't needed); and (3) the
CrautSolver
solves the system with a common matrix by the Craut method with partial
pivoting (the Craut method is the modification of Gauss eliminations). The
specialized solver, the
GreenSolver,
is used in spline approximation algorithm based on analytic representation of
spline. The
GeneralGreenSolver
extends the GreenSolver. It generalizes its superclass in solving of
SLAE with an arbitrary right-hand side.
RealCommonSolver
extends the RealDenseSolver. It is the basic class for solvers oriented
to the systems with common dense rectangular matrices. The system with square
matrix may be solved by a solver of this class also. This class supports the
balancing of matrices that may improve the accuracy of solution. The additional
solveT(y,x)
method allows to solve transposed system ATx = y using the factorization of
A matrix.
RotationSolver
based on Givens rotations; (2) the
ReflectionSolver
based on Hausholder reflections (it is theoretically 30% more efficient for
speed, but usually loses in accuracy; to provide the same accuracy as the
previous solver, the calculation of inner products with quadruple precision is
needed, but this is impossible in pure Java); and (3) the
EliminationSolver
based on Gauss eliminations with partial pivoting (it is twice faster that the
previous solver, but its accuracy is worst of all).
The ru.sscc.spline package contains basic spline approximation classes. We use the term spline for a function
|
|
A spline instance consists of 3 parts:
SplineBody interface.
SplineWorkspace is basic
class for a spline workspace.
The BaseSpline class is the basic class for all
splines. It implements common operations for all splines: the
dimension() operation returns the
dimension m of the space of points Rm for which the spline was
constructed; the
clone() operation creates a clone of the spline
(the clone shares the body and coefficients with the parent instance, but has
its own workspace); and the
flush() operation removes the workspace from
a spline. You can also use the
getBody() and
getContainer() operations for the
access to the spline body and coefficients.
All splines are serializable objects.
The BaseSpline instance cannot be created directly. But it has two subclasses having public constructors: the Spline and Splines. The first class performs operations with a unique spline, and the last class performs operations with many splines constructed on the same spline body.
The Spline class has a number of value()
operations. For 1D spline, you can use any of them, but, for multivariate
spline, you must use only the multi-D variant of the operation. If 1D spline
allows to calculate a spline derivative, the
value(point,derivative) operation
may be used. For the multivariate spline, you can calculate a spline value
only. If the calculation cannot be done for the parameters specified, the
IllegalArgumentException is thrown.
The Splines class has more value() operations: as
individual spline value calculation so as vector operations for all splines of
the instance.
The ru.sscc.spline.polynomial package is intended for the construction of 1D piecewise polynomial splines. The polynomial splines differ by the spline degree and spline representation. The construction of polynomial splines of odd degree with polynomial representation in mesh nodes is available now. We plan to add the construction of even degree splines and also implement the B-spline representation for all splines in future releases.
The POddSplineCreator class allows
to construct polynomial splines of odd degree. To describe its features, we
introduce the notation on the example of cubic splines.
Let x0 < x1 < ¼ < xn-1 be mesh nodes on the real axis and fi = f(xi) be the values of the function to be interpolated. The well-known cubic spline is the function constructed with the cubic polynomials defined at mesh cells (xi,xi+1), i = 0,¼,n-2, in such a way to provide its second derivative to be continuous function at mesh nodes. The interpolating cubic spline satisfies the interpolation conditions
|
|
To construct and use the natural interpolating cubic spline, write the
following code:
// Create a mesh array. It may be a double[], float[], or RealContainer instance.
// For example:
double[] mesh = new double[n];
// ... fill the mesh
// Create an array of interpolation values at mesh nodes. For example:
double[] values = new double[n];
// ... fill the array
// Create the natural interpolating cubic spline
Spline spl = POddSplineCreator.createSpline(2, mesh, values);
// Use the spline in calculations: x - a point on the real axis, f - the spline value
f = spl.value(x);
The first parameter of all createSpline methods in the
POddSplineCreator class
is the difference order that is tied with the spline degree by the
following rule: degree=2*order-1. So, to create the cubic spline, the
order must be equal to 2. If order=1, the piecewise linear spline will be
constructed. For order=3,4,..., the quintic spline, 7th degree spline, etc. will be constructed. The more order is used, the higher spline derivative is
continuous. Exactly, the (2*order-2)-th spline derivative is the
continuous function. The extrapolation (e.g. calculation of the spline
values for x < x0 and x > xn-1) is provided by the polynomials of
(order-1)-th degree. So, the cubic spline is extended by linear
functions, the quintic spline is extended by quadratic functions and so on.
The number of mesh nodes must satisfy the inequality n>=order. The mesh
nodes must be ordered in increasing or decreasing order and must be different.
If the mesh is incorrect, or some mesh nodes are very closed, or the difference
order is too large, the
CalculatingException will be thrown while
the spline construction. This means that the spline cannot be constructed. You
must catch this exception anywhere.
If the distance between mesh nodes is constant, you need not fill in the mesh
array. Use the following method instead:
Spline spl = POddSplineCreator.createSpline(order, x0, step, n, values);
Here x0 is the first mesh node, step is the distance between
neighbour mesh nodes, n is the number of mesh nodes, and values
is the array of interpolating values at mesh nodes.
When the function measurements are non-strict, the constructing of smoothing spline is more preferable than the interpolation. The smoothing spline goes "near" the function measurements. It is more smooth and may suppress unexpected overshoots in the measurements. The degree of smoothing is controlled by the smoothing parameter a > 0. When a tends to zero, the smoothing spline converges to the interpolation spline. When a tends to the positive infinity, the smoothing spline converges to the polynomial of the degree order-1 (this is the linear function for cubic spline). The problem is to select proper a providing the "best" smoothing results.
We select the smoothing parameter to satisfy the residual criterion
|
It isn't necessary to solve this equation with high accuracy. It is usually
enough to construct the spline, such that the residual belongs to the segment
[e(1-d),e(1+delta)]. Here d is the accuracy
level. The default accuracy level is set to 0.1. It may be changed by the
SmoothingPreparator.setAccuracy(delta)
method.
To construct a smoothing spline satisfying the residual criterion, the
following methods are used:
// For arbitrary mesh:
Spline spl = POddSplineCreator.createSpline(order, mesh, values, epsilon);
// For uniform mesh:
Spline spl = POddSplineCreator.createSpline(order, x0, step, n, values, epsilon);
If epsilon isn't less than emax, the spline will
coincide with the polynomial of the (order-1)-th degree-the solution
to the smoothing problem for a = ¥.
If you know the function measurement error for every mesh node, you'll want to construct the smoothing spline that may strongly deviate from the function measurements in "bad" nodes, and weekly deviate in the "good" nodes. To provide this, use the weighted residual criterion
|
To construct a smoothing spline satisfying the weighted residual criterion, the
following methods are used:
// For arbitrary mesh:
Spline spl = POddSplineCreator.createSpline(order, mesh, values, epsilon, weights);
// For uniform mesh:
Spline spl = POddSplineCreator.createSpline(order, x0, step, n, values, epsilon, weights);
Here the weights array contains positive weights assigned to mesh nodes.
If you want to construct many splines distinct in function measurements only, you must do the following:
POddSplineCreator class.
The Step 3 may be repeated many times.
Example 1 To construct a number of interpolating splines on an
arbitrary mesh, write down the following:
double[] mesh = new double[n];
// ... fill the mesh
OddSplineCreator creator = new POddSplineCreator(order, RealVector.by(mesh));
creator.prepareSolver(0); // Prepare solver for alpha=0
// Create array for the measurements
double[] values = new double[n];
// ... prepare function measurements in the array
// Construct the first spline:
Spline spl1 = creator.constructSpline(RealVector.by(values));
// ... prepare function measurements in the same array
// Construct the second spline:
Spline spl2 = creator.constructSpline(RealVector.by(values));
// And so on ...
The POddSplineCreator class is the implementation of the abstract
OddSplineCreator class. The conversion of the creator instance to its
superclass is preferable in this case, because the superclass has just the same
interface.
Example 2 To construct the weighted smoothing splines with the same
smoothing parameter a (1 in this example), write down the following:
OddSplineCreator creator = new POddSplineCreator(order, RealVector.by(mesh));
double[] weights = new double[n];
// ... fill the weight vector
creator.prepareWeights(weights);
creator.prepareSolver(1); // Prepare solver for alpha=1
// ... construct splines by the same manner as in the previous example
Example 3 To construct the smoothing spline with selection of
a from the residual criterion, you need not prepare the solver. Instead
of the preparation of the solver you can immediately use the constructing
method:
Spline spl1 = creator.constructSpline(RealVector.by(values), epsilon);
After that, the solver will be prepared. So, to construct other splines with
the same parameter a, you can call:
Spline spl2 = creator.constructSpline(RealVector.by(values));
// ...
Example 4 When many splines must be constructed in one Splines
instance, you should prepare the RealVectors instance containing all
functions measurements and use one of the following methods:
RealVectors measurements;
// ... prepare the measurements
// Construct splines by the rows of the measurements array:
Splines spls = creator.constructSplinesByRows(measurements);
// or construct splines by the columns of the measurements array:
Splines spls = creator.constructSplinesByColumns(measurements);
Example 5 If you want to create the Splines instance with
sequential calculation the splines' coefficients, do the following. Manually
construct the Splines instance:
Splines spls = new Splines(creator.getSplineBody(), count);
Here count is a number of splines to be calculated. After that, call the
calculate method in a cycle:
for(int i=0; i<count; i++) {
// ... prepare i-th measurements in the values array ...
// Calculate the coefficients of i-th spline
creator.calculate(RealVector.by(values), spls.getVector(i));
}
The ru.sscc.spline.analytic package is intended for the construction of multivariate analytic Duchon's splines. The Duchon's spline is the generalization of so-called 2D thin plate spline to multi-D case. It allows to construct spline for arbitrary number of independent variables on a scattered mesh of interpolating nodes. The Duchon's spline is constructed with the help of especial multivariate radial basic function g(x). The range of possible basic functions is wide, but we use 3 of them only:
The Duchon's spline has the form
|
The type of spline to be constructed is defined by the mode parameter. Four different modes may be used:
The more mode is used, the smoother spline is constructed.
The GSplineCreator class allows to
construct Duchon's splines. The use of this class is quite similar to the use
of the polynomial spline creator on a nonuniform mesh (see
Section 4.1) except two essential distinctions:
RealVectors or
ReducedMesh type.
The interpolation, fitting, and residual criterion fitting (weighted or not) are possible.
|
||||||||
| PREV NEXT | FRAMES NO FRAMES | |||||||