alglibinternal.h   alglibinternal.h 
skipping to change at line 115 skipping to change at line 115
///////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////
// //
// THIS SECTION CONTAINS COMPUTATIONAL CORE DECLARATIONS (FUNCTIONS) // THIS SECTION CONTAINS COMPUTATIONAL CORE DECLARATIONS (FUNCTIONS)
// //
///////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////
namespace alglib_impl namespace alglib_impl
{ {
ae_int_t getrdfserializationcode(ae_state *_state); ae_int_t getrdfserializationcode(ae_state *_state);
ae_int_t getkdtreeserializationcode(ae_state *_state); ae_int_t getkdtreeserializationcode(ae_state *_state);
ae_int_t getmlpserializationcode(ae_state *_state); ae_int_t getmlpserializationcode(ae_state *_state);
ae_int_t getmlpeserializationcode(ae_state *_state);
ae_int_t getrbfserializationcode(ae_state *_state);
ae_bool approxequalrel(double a, double b, double tol, ae_state *_state);
void taskgenint1d(double a, void taskgenint1d(double a,
double b, double b,
ae_int_t n, ae_int_t n,
/* Real */ ae_vector* x, /* Real */ ae_vector* x,
/* Real */ ae_vector* y, /* Real */ ae_vector* y,
ae_state *_state); ae_state *_state);
void taskgenint1dequidist(double a, void taskgenint1dequidist(double a,
double b, double b,
ae_int_t n, ae_int_t n,
/* Real */ ae_vector* x, /* Real */ ae_vector* x,
skipping to change at line 142 skipping to change at line 145
ae_state *_state); ae_state *_state);
void taskgenint1dcheb2(double a, void taskgenint1dcheb2(double a,
double b, double b,
ae_int_t n, ae_int_t n,
/* Real */ ae_vector* x, /* Real */ ae_vector* x,
/* Real */ ae_vector* y, /* Real */ ae_vector* y,
ae_state *_state); ae_state *_state);
ae_bool aredistinct(/* Real */ ae_vector* x, ae_bool aredistinct(/* Real */ ae_vector* x,
ae_int_t n, ae_int_t n,
ae_state *_state); ae_state *_state);
ae_bool aresameboolean(ae_bool v1, ae_bool v2, ae_state *_state);
void bvectorsetlengthatleast(/* Boolean */ ae_vector* x, void bvectorsetlengthatleast(/* Boolean */ ae_vector* x,
ae_int_t n, ae_int_t n,
ae_state *_state); ae_state *_state);
void ivectorsetlengthatleast(/* Integer */ ae_vector* x, void ivectorsetlengthatleast(/* Integer */ ae_vector* x,
ae_int_t n, ae_int_t n,
ae_state *_state); ae_state *_state);
void rvectorsetlengthatleast(/* Real */ ae_vector* x, void rvectorsetlengthatleast(/* Real */ ae_vector* x,
ae_int_t n, ae_int_t n,
ae_state *_state); ae_state *_state);
void rmatrixsetlengthatleast(/* Real */ ae_matrix* x, void rmatrixsetlengthatleast(/* Real */ ae_matrix* x,
skipping to change at line 279 skipping to change at line 283
void tagsortfastr(/* Real */ ae_vector* a, void tagsortfastr(/* Real */ ae_vector* a,
/* Real */ ae_vector* b, /* Real */ ae_vector* b,
/* Real */ ae_vector* bufa, /* Real */ ae_vector* bufa,
/* Real */ ae_vector* bufb, /* Real */ ae_vector* bufb,
ae_int_t n, ae_int_t n,
ae_state *_state); ae_state *_state);
void tagsortfast(/* Real */ ae_vector* a, void tagsortfast(/* Real */ ae_vector* a,
/* Real */ ae_vector* bufa, /* Real */ ae_vector* bufa,
ae_int_t n, ae_int_t n,
ae_state *_state); ae_state *_state);
void tagsortmiddleir(/* Integer */ ae_vector* a,
/* Real */ ae_vector* b,
ae_int_t offset,
ae_int_t n,
ae_state *_state);
void tagheappushi(/* Real */ ae_vector* a, void tagheappushi(/* Real */ ae_vector* a,
/* Integer */ ae_vector* b, /* Integer */ ae_vector* b,
ae_int_t* n, ae_int_t* n,
double va, double va,
ae_int_t vb, ae_int_t vb,
ae_state *_state); ae_state *_state);
void tagheapreplacetopi(/* Real */ ae_vector* a, void tagheapreplacetopi(/* Real */ ae_vector* a,
/* Integer */ ae_vector* b, /* Integer */ ae_vector* b,
ae_int_t n, ae_int_t n,
double va, double va,
ae_int_t vb, ae_int_t vb,
ae_state *_state); ae_state *_state);
void tagheappopi(/* Real */ ae_vector* a, void tagheappopi(/* Real */ ae_vector* a,
/* Integer */ ae_vector* b, /* Integer */ ae_vector* b,
ae_int_t* n, ae_int_t* n,
ae_state *_state); ae_state *_state);
ae_int_t lowerbound(/* Real */ ae_vector* a,
ae_int_t n,
double t,
ae_state *_state);
ae_int_t upperbound(/* Real */ ae_vector* a,
ae_int_t n,
double t,
ae_state *_state);
void rankx(/* Real */ ae_vector* x, void rankx(/* Real */ ae_vector* x,
ae_int_t n, ae_int_t n,
apbuffers* buf, apbuffers* buf,
ae_state *_state); ae_state *_state);
ae_bool cmatrixrank1f(ae_int_t m, ae_bool cmatrixrank1f(ae_int_t m,
ae_int_t n, ae_int_t n,
/* Complex */ ae_matrix* a, /* Complex */ ae_matrix* a,
ae_int_t ia, ae_int_t ia,
ae_int_t ja, ae_int_t ja,
/* Complex */ ae_vector* u, /* Complex */ ae_vector* u,
skipping to change at line 449 skipping to change at line 466
ae_int_t optypea, ae_int_t optypea,
/* Complex */ ae_matrix* b, /* Complex */ ae_matrix* b,
ae_int_t ib, ae_int_t ib,
ae_int_t jb, ae_int_t jb,
ae_int_t optypeb, ae_int_t optypeb,
ae_complex beta, ae_complex beta,
/* Complex */ ae_matrix* c, /* Complex */ ae_matrix* c,
ae_int_t ic, ae_int_t ic,
ae_int_t jc, ae_int_t jc,
ae_state *_state); ae_state *_state);
void hermitianmatrixvectormultiply(/* Complex */ ae_matrix* a,
ae_bool isupper,
ae_int_t i1,
ae_int_t i2,
/* Complex */ ae_vector* x,
ae_complex alpha,
/* Complex */ ae_vector* y,
ae_state *_state);
void hermitianrank2update(/* Complex */ ae_matrix* a,
ae_bool isupper,
ae_int_t i1,
ae_int_t i2,
/* Complex */ ae_vector* x,
/* Complex */ ae_vector* y,
/* Complex */ ae_vector* t,
ae_complex alpha,
ae_state *_state);
void generatereflection(/* Real */ ae_vector* x,
ae_int_t n,
double* tau,
ae_state *_state);
void applyreflectionfromtheleft(/* Real */ ae_matrix* c,
double tau,
/* Real */ ae_vector* v,
ae_int_t m1,
ae_int_t m2,
ae_int_t n1,
ae_int_t n2,
/* Real */ ae_vector* work,
ae_state *_state);
void applyreflectionfromtheright(/* Real */ ae_matrix* c,
double tau,
/* Real */ ae_vector* v,
ae_int_t m1,
ae_int_t m2,
ae_int_t n1,
ae_int_t n2,
/* Real */ ae_vector* work,
ae_state *_state);
void complexgeneratereflection(/* Complex */ ae_vector* x,
ae_int_t n,
ae_complex* tau,
ae_state *_state);
void complexapplyreflectionfromtheleft(/* Complex */ ae_matrix* c,
ae_complex tau,
/* Complex */ ae_vector* v,
ae_int_t m1,
ae_int_t m2,
ae_int_t n1,
ae_int_t n2,
/* Complex */ ae_vector* work,
ae_state *_state);
void complexapplyreflectionfromtheright(/* Complex */ ae_matrix* c,
ae_complex tau,
/* Complex */ ae_vector* v,
ae_int_t m1,
ae_int_t m2,
ae_int_t n1,
ae_int_t n2,
/* Complex */ ae_vector* work,
ae_state *_state);
void symmetricmatrixvectormultiply(/* Real */ ae_matrix* a,
ae_bool isupper,
ae_int_t i1,
ae_int_t i2,
/* Real */ ae_vector* x,
double alpha,
/* Real */ ae_vector* y,
ae_state *_state);
void symmetricrank2update(/* Real */ ae_matrix* a,
ae_bool isupper,
ae_int_t i1,
ae_int_t i2,
/* Real */ ae_vector* x,
/* Real */ ae_vector* y,
/* Real */ ae_vector* t,
double alpha,
ae_state *_state);
double vectornorm2(/* Real */ ae_vector* x, double vectornorm2(/* Real */ ae_vector* x,
ae_int_t i1, ae_int_t i1,
ae_int_t i2, ae_int_t i2,
ae_state *_state); ae_state *_state);
ae_int_t vectoridxabsmax(/* Real */ ae_vector* x, ae_int_t vectoridxabsmax(/* Real */ ae_vector* x,
ae_int_t i1, ae_int_t i1,
ae_int_t i2, ae_int_t i2,
ae_state *_state); ae_state *_state);
ae_int_t columnidxabsmax(/* Real */ ae_matrix* x, ae_int_t columnidxabsmax(/* Real */ ae_matrix* x,
ae_int_t i1, ae_int_t i1,
skipping to change at line 540 skipping to change at line 635
ae_bool transb, ae_bool transb,
double alpha, double alpha,
/* Real */ ae_matrix* c, /* Real */ ae_matrix* c,
ae_int_t ci1, ae_int_t ci1,
ae_int_t ci2, ae_int_t ci2,
ae_int_t cj1, ae_int_t cj1,
ae_int_t cj2, ae_int_t cj2,
double beta, double beta,
/* Real */ ae_vector* work, /* Real */ ae_vector* work,
ae_state *_state); ae_state *_state);
void hermitianmatrixvectormultiply(/* Complex */ ae_matrix* a,
ae_bool isupper,
ae_int_t i1,
ae_int_t i2,
/* Complex */ ae_vector* x,
ae_complex alpha,
/* Complex */ ae_vector* y,
ae_state *_state);
void hermitianrank2update(/* Complex */ ae_matrix* a,
ae_bool isupper,
ae_int_t i1,
ae_int_t i2,
/* Complex */ ae_vector* x,
/* Complex */ ae_vector* y,
/* Complex */ ae_vector* t,
ae_complex alpha,
ae_state *_state);
void generatereflection(/* Real */ ae_vector* x,
ae_int_t n,
double* tau,
ae_state *_state);
void applyreflectionfromtheleft(/* Real */ ae_matrix* c,
double tau,
/* Real */ ae_vector* v,
ae_int_t m1,
ae_int_t m2,
ae_int_t n1,
ae_int_t n2,
/* Real */ ae_vector* work,
ae_state *_state);
void applyreflectionfromtheright(/* Real */ ae_matrix* c,
double tau,
/* Real */ ae_vector* v,
ae_int_t m1,
ae_int_t m2,
ae_int_t n1,
ae_int_t n2,
/* Real */ ae_vector* work,
ae_state *_state);
void complexgeneratereflection(/* Complex */ ae_vector* x,
ae_int_t n,
ae_complex* tau,
ae_state *_state);
void complexapplyreflectionfromtheleft(/* Complex */ ae_matrix* c,
ae_complex tau,
/* Complex */ ae_vector* v,
ae_int_t m1,
ae_int_t m2,
ae_int_t n1,
ae_int_t n2,
/* Complex */ ae_vector* work,
ae_state *_state);
void complexapplyreflectionfromtheright(/* Complex */ ae_matrix* c,
ae_complex tau,
/* Complex */ ae_vector* v,
ae_int_t m1,
ae_int_t m2,
ae_int_t n1,
ae_int_t n2,
/* Complex */ ae_vector* work,
ae_state *_state);
void symmetricmatrixvectormultiply(/* Real */ ae_matrix* a,
ae_bool isupper,
ae_int_t i1,
ae_int_t i2,
/* Real */ ae_vector* x,
double alpha,
/* Real */ ae_vector* y,
ae_state *_state);
void symmetricrank2update(/* Real */ ae_matrix* a,
ae_bool isupper,
ae_int_t i1,
ae_int_t i2,
/* Real */ ae_vector* x,
/* Real */ ae_vector* y,
/* Real */ ae_vector* t,
double alpha,
ae_state *_state);
void applyrotationsfromtheleft(ae_bool isforward, void applyrotationsfromtheleft(ae_bool isforward,
ae_int_t m1, ae_int_t m1,
ae_int_t m2, ae_int_t m2,
ae_int_t n1, ae_int_t n1,
ae_int_t n2, ae_int_t n2,
/* Real */ ae_vector* c, /* Real */ ae_vector* c,
/* Real */ ae_vector* s, /* Real */ ae_vector* s,
/* Real */ ae_matrix* a, /* Real */ ae_matrix* a,
/* Real */ ae_vector* work, /* Real */ ae_vector* work,
ae_state *_state); ae_state *_state);
skipping to change at line 746 skipping to change at line 763
ae_int_t* info, ae_int_t* info,
double* stp, double* stp,
double* f, double* f,
ae_state *_state); ae_state *_state);
ae_bool _linminstate_init(linminstate* p, ae_state *_state, ae_bool make_au tomatic); ae_bool _linminstate_init(linminstate* p, ae_state *_state, ae_bool make_au tomatic);
ae_bool _linminstate_init_copy(linminstate* dst, linminstate* src, ae_state *_state, ae_bool make_automatic); ae_bool _linminstate_init_copy(linminstate* dst, linminstate* src, ae_state *_state, ae_bool make_automatic);
void _linminstate_clear(linminstate* p); void _linminstate_clear(linminstate* p);
ae_bool _armijostate_init(armijostate* p, ae_state *_state, ae_bool make_au tomatic); ae_bool _armijostate_init(armijostate* p, ae_state *_state, ae_bool make_au tomatic);
ae_bool _armijostate_init_copy(armijostate* dst, armijostate* src, ae_state *_state, ae_bool make_automatic); ae_bool _armijostate_init_copy(armijostate* dst, armijostate* src, ae_state *_state, ae_bool make_automatic);
void _armijostate_clear(armijostate* p); void _armijostate_clear(armijostate* p);
void trimprepare(double f, double* threshold, ae_state *_state);
void trimfunction(double* f,
/* Real */ ae_vector* g,
ae_int_t n,
double threshold,
ae_state *_state);
void ftbasegeneratecomplexfftplan(ae_int_t n, void ftbasegeneratecomplexfftplan(ae_int_t n,
ftplan* plan, ftplan* plan,
ae_state *_state); ae_state *_state);
void ftbasegeneraterealfftplan(ae_int_t n, ftplan* plan, ae_state *_state); void ftbasegeneraterealfftplan(ae_int_t n, ftplan* plan, ae_state *_state);
void ftbasegeneraterealfhtplan(ae_int_t n, ftplan* plan, ae_state *_state); void ftbasegeneraterealfhtplan(ae_int_t n, ftplan* plan, ae_state *_state);
void ftbaseexecuteplan(/* Real */ ae_vector* a, void ftbaseexecuteplan(/* Real */ ae_vector* a,
ae_int_t aoffset, ae_int_t aoffset,
ae_int_t n, ae_int_t n,
ftplan* plan, ftplan* plan,
ae_state *_state); ae_state *_state);
 End of changes. 7 change blocks. 
84 lines changed or deleted 95 lines changed or added


 alglibmisc.h   alglibmisc.h 
skipping to change at line 218 skipping to change at line 218
Random number generator: exponential distribution Random number generator: exponential distribution
State structure must be initialized with HQRNDRandomize() or HQRNDSeed(). State structure must be initialized with HQRNDRandomize() or HQRNDSeed().
-- ALGLIB -- -- ALGLIB --
Copyright 11.08.2007 by Bochkanov Sergey Copyright 11.08.2007 by Bochkanov Sergey
*************************************************************************/ *************************************************************************/
double hqrndexponential(const hqrndstate &state, const double lambdav); double hqrndexponential(const hqrndstate &state, const double lambdav);
/************************************************************************* /*************************************************************************
This function generates random number from discrete distribution given by
finite sample X.
INPUT PARAMETERS
State - high quality random number generator, must be
initialized with HQRNDRandomize() or HQRNDSeed().
X - finite sample
N - number of elements to use, N>=1
RESULT
this function returns one of the X[i] for random i=0..N-1
-- ALGLIB --
Copyright 08.11.2011 by Bochkanov Sergey
*************************************************************************/
double hqrnddiscrete(const hqrndstate &state, const real_1d_array &x, const
ae_int_t n);
/*************************************************************************
This function generates random number from continuous distribution given
by finite sample X.
INPUT PARAMETERS
State - high quality random number generator, must be
initialized with HQRNDRandomize() or HQRNDSeed().
X - finite sample, array[N] (can be larger, in this case only
leading N elements are used). THIS ARRAY MUST BE SORTED BY
ASCENDING.
N - number of elements to use, N>=1
RESULT
this function returns random number from continuous distribution which
tries to approximate X as mush as possible. min(X)<=Result<=max(X).
-- ALGLIB --
Copyright 08.11.2011 by Bochkanov Sergey
*************************************************************************/
double hqrndcontinuous(const hqrndstate &state, const real_1d_array &x, con
st ae_int_t n);
/*************************************************************************
This function serializes data structure to string. This function serializes data structure to string.
Important properties of s_out: Important properties of s_out:
* it contains alphanumeric characters, dots, underscores, minus signs * it contains alphanumeric characters, dots, underscores, minus signs
* these symbols are grouped into words, which are separated by spaces * these symbols are grouped into words, which are separated by spaces
and Windows-style (CR+LF) newlines and Windows-style (CR+LF) newlines
* although serializer uses spaces and CR+LF as separators, you can * although serializer uses spaces and CR+LF as separators, you can
replace any separator character by arbitrary combination of spaces, replace any separator character by arbitrary combination of spaces,
tabs, Windows or Unix newlines. It allows flexible reformatting of tabs, Windows or Unix newlines. It allows flexible reformatting of
the string in case you want to include it into text or XML file. the string in case you want to include it into text or XML file.
skipping to change at line 254 skipping to change at line 293
/************************************************************************* /*************************************************************************
KD-tree creation KD-tree creation
This subroutine creates KD-tree from set of X-values and optional Y-values This subroutine creates KD-tree from set of X-values and optional Y-values
INPUT PARAMETERS INPUT PARAMETERS
XY - dataset, array[0..N-1,0..NX+NY-1]. XY - dataset, array[0..N-1,0..NX+NY-1].
one row corresponds to one point. one row corresponds to one point.
first NX columns contain X-values, next NY (NY may be zero) first NX columns contain X-values, next NY (NY may be zero)
columns may contain associated Y-values columns may contain associated Y-values
N - number of points, N>=1 N - number of points, N>=0.
NX - space dimension, NX>=1. NX - space dimension, NX>=1.
NY - number of optional Y-values, NY>=0. NY - number of optional Y-values, NY>=0.
NormType- norm type: NormType- norm type:
* 0 denotes infinity-norm * 0 denotes infinity-norm
* 1 denotes 1-norm * 1 denotes 1-norm
* 2 denotes 2-norm (Euclidean norm) * 2 denotes 2-norm (Euclidean norm)
OUTPUT PARAMETERS OUTPUT PARAMETERS
KDT - KD-tree KDT - KD-tree
skipping to change at line 294 skipping to change at line 333
This subroutine creates KD-tree from set of X-values, integer tags and This subroutine creates KD-tree from set of X-values, integer tags and
optional Y-values optional Y-values
INPUT PARAMETERS INPUT PARAMETERS
XY - dataset, array[0..N-1,0..NX+NY-1]. XY - dataset, array[0..N-1,0..NX+NY-1].
one row corresponds to one point. one row corresponds to one point.
first NX columns contain X-values, next NY (NY may be zero) first NX columns contain X-values, next NY (NY may be zero)
columns may contain associated Y-values columns may contain associated Y-values
Tags - tags, array[0..N-1], contains integer tags associated Tags - tags, array[0..N-1], contains integer tags associated
with points. with points.
N - number of points, N>=1 N - number of points, N>=0
NX - space dimension, NX>=1. NX - space dimension, NX>=1.
NY - number of optional Y-values, NY>=0. NY - number of optional Y-values, NY>=0.
NormType- norm type: NormType- norm type:
* 0 denotes infinity-norm * 0 denotes infinity-norm
* 1 denotes 1-norm * 1 denotes 1-norm
* 2 denotes 2-norm (Euclidean norm) * 2 denotes 2-norm (Euclidean norm)
OUTPUT PARAMETERS OUTPUT PARAMETERS
KDT - KD-tree KDT - KD-tree
skipping to change at line 627 skipping to change at line 666
ae_int_t hqrnduniformi(hqrndstate* state, ae_int_t n, ae_state *_state); ae_int_t hqrnduniformi(hqrndstate* state, ae_int_t n, ae_state *_state);
double hqrndnormal(hqrndstate* state, ae_state *_state); double hqrndnormal(hqrndstate* state, ae_state *_state);
void hqrndunit2(hqrndstate* state, double* x, double* y, ae_state *_state); void hqrndunit2(hqrndstate* state, double* x, double* y, ae_state *_state);
void hqrndnormal2(hqrndstate* state, void hqrndnormal2(hqrndstate* state,
double* x1, double* x1,
double* x2, double* x2,
ae_state *_state); ae_state *_state);
double hqrndexponential(hqrndstate* state, double hqrndexponential(hqrndstate* state,
double lambdav, double lambdav,
ae_state *_state); ae_state *_state);
double hqrnddiscrete(hqrndstate* state,
/* Real */ ae_vector* x,
ae_int_t n,
ae_state *_state);
double hqrndcontinuous(hqrndstate* state,
/* Real */ ae_vector* x,
ae_int_t n,
ae_state *_state);
ae_bool _hqrndstate_init(hqrndstate* p, ae_state *_state, ae_bool make_auto matic); ae_bool _hqrndstate_init(hqrndstate* p, ae_state *_state, ae_bool make_auto matic);
ae_bool _hqrndstate_init_copy(hqrndstate* dst, hqrndstate* src, ae_state *_ state, ae_bool make_automatic); ae_bool _hqrndstate_init_copy(hqrndstate* dst, hqrndstate* src, ae_state *_ state, ae_bool make_automatic);
void _hqrndstate_clear(hqrndstate* p); void _hqrndstate_clear(hqrndstate* p);
void kdtreebuild(/* Real */ ae_matrix* xy, void kdtreebuild(/* Real */ ae_matrix* xy,
ae_int_t n, ae_int_t n,
ae_int_t nx, ae_int_t nx,
ae_int_t ny, ae_int_t ny,
ae_int_t normtype, ae_int_t normtype,
kdtree* kdt, kdtree* kdt,
ae_state *_state); ae_state *_state);
 End of changes. 4 change blocks. 
2 lines changed or deleted 51 lines changed or added


 ap.h   ap.h 
skipping to change at line 532 skipping to change at line 532
double ae_sin(double x, ae_state *state); double ae_sin(double x, ae_state *state);
double ae_cos(double x, ae_state *state); double ae_cos(double x, ae_state *state);
double ae_tan(double x, ae_state *state); double ae_tan(double x, ae_state *state);
double ae_sinh(double x, ae_state *state); double ae_sinh(double x, ae_state *state);
double ae_cosh(double x, ae_state *state); double ae_cosh(double x, ae_state *state);
double ae_tanh(double x, ae_state *state); double ae_tanh(double x, ae_state *state);
double ae_asin(double x, ae_state *state); double ae_asin(double x, ae_state *state);
double ae_acos(double x, ae_state *state); double ae_acos(double x, ae_state *state);
double ae_atan(double x, ae_state *state); double ae_atan(double x, ae_state *state);
double ae_atan2(double x, double y, ae_state *state); double ae_atan2(double y, double x, ae_state *state);
double ae_log(double x, ae_state *state); double ae_log(double x, ae_state *state);
double ae_pow(double x, double y, ae_state *state); double ae_pow(double x, double y, ae_state *state);
double ae_exp(double x, ae_state *state); double ae_exp(double x, ae_state *state);
/************************************************************************ /************************************************************************
Complex math functions: Complex math functions:
* basic arithmetic operations * basic arithmetic operations
* standard functions * standard functions
************************************************************************/ ************************************************************************/
skipping to change at line 631 skipping to change at line 631
} rcommstate; } rcommstate;
ae_bool _rcommstate_init(rcommstate* p, ae_state *_state, ae_bool make_auto matic); ae_bool _rcommstate_init(rcommstate* p, ae_state *_state, ae_bool make_auto matic);
ae_bool _rcommstate_init_copy(rcommstate* dst, rcommstate* src, ae_state *_ state, ae_bool make_automatic); ae_bool _rcommstate_init_copy(rcommstate* dst, rcommstate* src, ae_state *_ state, ae_bool make_automatic);
void _rcommstate_clear(rcommstate* p); void _rcommstate_clear(rcommstate* p);
#ifdef AE_USE_ALLOC_COUNTER #ifdef AE_USE_ALLOC_COUNTER
extern ae_int64_t _alloc_counter; extern ae_int64_t _alloc_counter;
#endif #endif
/************************************************************************ /************************************************************************
debug functions (must be turned on by preprocessor definitions: debug functions (must be turned on by preprocessor definitions):
* tickcount(), which is wrapper around GetTickCount() * tickcount(), which is wrapper around GetTickCount()
* flushconsole(), fluches console
* ae_debugrng(), returns random number generated with high-quality random n
umbers generator
* ae_set_seed(), sets seed of the debug RNG (NON-THREAD-SAFE!!!)
* ae_get_seed(), returns two seed values of the debug RNG (NON-THREAD-SAFE!
!!)
************************************************************************/ ************************************************************************/
#ifdef AE_DEBUG4WINDOWS #ifdef AE_DEBUG4WINDOWS
#include <windows.h> #include <windows.h>
#include <stdio.h> #include <stdio.h>
#define tickcount(s) GetTickCount() #define tickcount(s) GetTickCount()
#define flushconsole(s) fflush(stdout) #define flushconsole(s) fflush(stdout)
#endif #endif
#ifdef AE_DEBUGRNG
ae_int_t ae_debugrng();
void ae_set_seed(ae_int_t s0, ae_int_t s1);
void ae_get_seed(ae_int_t *s0, ae_int_t *s1);
#endif
} }
///////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////
// //
// THIS SECTION CONTAINS DECLARATIONS FOR C++ RELATED FUNCTIONALITY // THIS SECTION CONTAINS DECLARATIONS FOR C++ RELATED FUNCTIONALITY
// //
///////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////
namespace alglib namespace alglib
skipping to change at line 1069 skipping to change at line 1078
extern const double machineepsilon; extern const double machineepsilon;
extern const double maxrealnumber; extern const double maxrealnumber;
extern const double minrealnumber; extern const double minrealnumber;
extern const double fp_nan; extern const double fp_nan;
extern const double fp_posinf; extern const double fp_posinf;
extern const double fp_neginf; extern const double fp_neginf;
extern const ae_int_t endianness; extern const ae_int_t endianness;
int sign(double x); int sign(double x);
double randomreal(); double randomreal();
int randominteger(int maxv); ae_int_t randominteger(ae_int_t maxv);
int round(double x); int round(double x);
int trunc(double x); int trunc(double x);
int ifloor(double x); int ifloor(double x);
int iceil(double x); int iceil(double x);
double pi(); double pi();
double sqr(double x); double sqr(double x);
int maxint(int m1, int m2); int maxint(int m1, int m2);
int minint(int m1, int m2); int minint(int m1, int m2);
double maxreal(double m1, double m2); double maxreal(double m1, double m2);
double minreal(double m1, double m2); double minreal(double m1, double m2);
 End of changes. 5 change blocks. 
3 lines changed or deleted 14 lines changed or added


 dataanalysis.h   dataanalysis.h 
skipping to change at line 25 skipping to change at line 25
A copy of the GNU General Public License is available at A copy of the GNU General Public License is available at
http://www.fsf.org/licensing/licenses http://www.fsf.org/licensing/licenses
>>> END OF LICENSE >>> >>> END OF LICENSE >>>
*************************************************************************/ *************************************************************************/
#ifndef _dataanalysis_pkg_h #ifndef _dataanalysis_pkg_h
#define _dataanalysis_pkg_h #define _dataanalysis_pkg_h
#include "ap.h" #include "ap.h"
#include "alglibinternal.h" #include "alglibinternal.h"
#include "linalg.h" #include "linalg.h"
#include "statistics.h" #include "statistics.h"
#include "alglibmisc.h"
#include "specialfunctions.h" #include "specialfunctions.h"
#include "alglibmisc.h"
#include "solvers.h" #include "solvers.h"
#include "optimization.h" #include "optimization.h"
///////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////
// //
// THIS SECTION CONTAINS COMPUTATIONAL CORE DECLARATIONS (DATATYPES) // THIS SECTION CONTAINS COMPUTATIONAL CORE DECLARATIONS (DATATYPES)
// //
///////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////
namespace alglib_impl namespace alglib_impl
{ {
skipping to change at line 207 skipping to change at line 207
typedef struct typedef struct
{ {
double relclserror; double relclserror;
double avgce; double avgce;
double rmserror; double rmserror;
double avgerror; double avgerror;
double avgrelerror; double avgrelerror;
} mlpcvreport; } mlpcvreport;
typedef struct typedef struct
{ {
ae_vector structinfo;
ae_int_t ensemblesize; ae_int_t ensemblesize;
ae_int_t nin;
ae_int_t nout;
ae_int_t wcount;
ae_bool issoftmax;
ae_bool postprocessing;
ae_vector weights; ae_vector weights;
ae_vector columnmeans; ae_vector columnmeans;
ae_vector columnsigmas; ae_vector columnsigmas;
ae_int_t serializedlen; multilayerperceptron network;
ae_vector serializedmlp;
ae_vector tmpweights;
ae_vector tmpmeans;
ae_vector tmpsigmas;
ae_vector neurons;
ae_vector dfdnet;
ae_vector y; ae_vector y;
} mlpensemble; } mlpensemble;
} }
///////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////
// //
// THIS SECTION CONTAINS C++ INTERFACE // THIS SECTION CONTAINS C++ INTERFACE
// //
///////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////
skipping to change at line 882 skipping to change at line 870
Its meaning for regression task is obvious. As for Its meaning for regression task is obvious. As for
classification task, it means average relative error when estimating classification task, it means average relative error when estimating
posterior probability of belonging to the correct class. posterior probability of belonging to the correct class.
-- ALGLIB -- -- ALGLIB --
Copyright 16.02.2009 by Bochkanov Sergey Copyright 16.02.2009 by Bochkanov Sergey
*************************************************************************/ *************************************************************************/
double dfavgrelerror(const decisionforest &df, const real_2d_array &xy, con st ae_int_t npoints); double dfavgrelerror(const decisionforest &df, const real_2d_array &xy, con st ae_int_t npoints);
/************************************************************************* /*************************************************************************
k-means++ clusterization
INPUT PARAMETERS:
XY - dataset, array [0..NPoints-1,0..NVars-1].
NPoints - dataset size, NPoints>=K
NVars - number of variables, NVars>=1
K - desired number of clusters, K>=1
Restarts - number of restarts, Restarts>=1
OUTPUT PARAMETERS:
Info - return code:
* -3, if task is degenerate (number of distinct points
is
less than K)
* -1, if incorrect NPoints/NFeatures/K/Restarts was pas
sed
* 1, if subroutine finished successfully
C - array[0..NVars-1,0..K-1].matrix whose columns store
cluster's centers
XYC - array[NPoints], which contains cluster indexes
-- ALGLIB --
Copyright 21.03.2009 by Bochkanov Sergey
*************************************************************************/
void kmeansgenerate(const real_2d_array &xy, const ae_int_t npoints, const
ae_int_t nvars, const ae_int_t k, const ae_int_t restarts, ae_int_t &info,
real_2d_array &c, integer_1d_array &xyc);
/*************************************************************************
Multiclass Fisher LDA
Subroutine finds coefficients of linear combination which optimally separat
es
training set on classes.
INPUT PARAMETERS:
XY - training set, array[0..NPoints-1,0..NVars].
First NVars columns store values of independent
variables, next column stores number of class (from 0
to NClasses-1) which dataset element belongs to. Fracti
onal
values are rounded to nearest integer.
NPoints - training set size, NPoints>=0
NVars - number of independent variables, NVars>=1
NClasses - number of classes, NClasses>=2
OUTPUT PARAMETERS:
Info - return code:
* -4, if internal EVD subroutine hasn't converged
* -2, if there is a point with class number
outside of [0..NClasses-1].
* -1, if incorrect parameters was passed (NPoints<0,
NVars<1, NClasses<2)
* 1, if task has been solved
* 2, if there was a multicollinearity in training set,
but task has been solved.
W - linear combination coefficients, array[0..NVars-1]
-- ALGLIB --
Copyright 31.05.2008 by Bochkanov Sergey
*************************************************************************/
void fisherlda(const real_2d_array &xy, const ae_int_t npoints, const ae_in
t_t nvars, const ae_int_t nclasses, ae_int_t &info, real_1d_array &w);
/*************************************************************************
N-dimensional multiclass Fisher LDA
Subroutine finds coefficients of linear combinations which optimally separa
tes
training set on classes. It returns N-dimensional basis whose vector are so
rted
by quality of training set separation (in descending order).
INPUT PARAMETERS:
XY - training set, array[0..NPoints-1,0..NVars].
First NVars columns store values of independent
variables, next column stores number of class (from 0
to NClasses-1) which dataset element belongs to. Fracti
onal
values are rounded to nearest integer.
NPoints - training set size, NPoints>=0
NVars - number of independent variables, NVars>=1
NClasses - number of classes, NClasses>=2
OUTPUT PARAMETERS:
Info - return code:
* -4, if internal EVD subroutine hasn't converged
* -2, if there is a point with class number
outside of [0..NClasses-1].
* -1, if incorrect parameters was passed (NPoints<0,
NVars<1, NClasses<2)
* 1, if task has been solved
* 2, if there was a multicollinearity in training set,
but task has been solved.
W - basis, array[0..NVars-1,0..NVars-1]
columns of matrix stores basis vectors, sorted by
quality of training set separation (in descending order
)
-- ALGLIB --
Copyright 31.05.2008 by Bochkanov Sergey
*************************************************************************/
void fisherldan(const real_2d_array &xy, const ae_int_t npoints, const ae_i
nt_t nvars, const ae_int_t nclasses, ae_int_t &info, real_2d_array &w);
/*************************************************************************
Linear regression Linear regression
Subroutine builds model: Subroutine builds model:
Y = A(0)*X[0] + ... + A(N-1)*X[N-1] + A(N) Y = A(0)*X[0] + ... + A(N-1)*X[N-1] + A(N)
and model found in ALGLIB format, covariation matrix, training set errors and model found in ALGLIB format, covariation matrix, training set errors
(rms, average, average relative) and leave-one-out cross-validation (rms, average, average relative) and leave-one-out cross-validation
estimate of the generalization error. CV estimate calculated using fast estimate of the generalization error. CV estimate calculated using fast
algorithm with O(NPoints*NVars) complexity. algorithm with O(NPoints*NVars) complexity.
skipping to change at line 1164 skipping to change at line 1058
RESULT: RESULT:
average relative error. average relative error.
-- ALGLIB -- -- ALGLIB --
Copyright 30.08.2008 by Bochkanov Sergey Copyright 30.08.2008 by Bochkanov Sergey
*************************************************************************/ *************************************************************************/
double lravgrelerror(const linearmodel &lm, const real_2d_array &xy, const ae_int_t npoints); double lravgrelerror(const linearmodel &lm, const real_2d_array &xy, const ae_int_t npoints);
/************************************************************************* /*************************************************************************
Filters: simple moving averages (unsymmetric).
This filter replaces array by results of SMA(K) filter. SMA(K) is defined
as filter which averages at most K previous points (previous - not points
AROUND central point) - or less, in case of the first K-1 points.
INPUT PARAMETERS:
X - array[N], array to process. It can be larger than N,
in this case only first N points are processed.
N - points count, N>=0
K - K>=1 (K can be larger than N , such cases will be
correctly handled). Window width. K=1 corresponds to
identity transformation (nothing changes).
OUTPUT PARAMETERS:
X - array, whose first N elements were processed with SMA(K
)
NOTE 1: this function uses efficient in-place algorithm which does not
allocate temporary arrays.
NOTE 2: this algorithm makes only one pass through array and uses running
sum to speed-up calculation of the averages. Additional measures
are taken to ensure that running sum on a long sequence of zero
elements will be correctly reset to zero even in the presence of
round-off error.
NOTE 3: this is unsymmetric version of the algorithm, which does NOT
averages points after the current one. Only X[i], X[i-1], ... are
used when calculating new value of X[i]. We should also note that
this algorithm uses BOTH previous points and current one, i.e.
new value of X[i] depends on BOTH previous point and X[i] itself.
-- ALGLIB --
Copyright 25.10.2011 by Bochkanov Sergey
*************************************************************************/
void filtersma(real_1d_array &x, const ae_int_t n, const ae_int_t k);
void filtersma(real_1d_array &x, const ae_int_t k);
/*************************************************************************
Filters: exponential moving averages.
This filter replaces array by results of EMA(alpha) filter. EMA(alpha) is
defined as filter which replaces X[] by S[]:
S[0] = X[0]
S[t] = alpha*X[t] + (1-alpha)*S[t-1]
INPUT PARAMETERS:
X - array[N], array to process. It can be larger than N,
in this case only first N points are processed.
N - points count, N>=0
alpha - 0<alpha<=1, smoothing parameter.
OUTPUT PARAMETERS:
X - array, whose first N elements were processed
with EMA(alpha)
NOTE 1: this function uses efficient in-place algorithm which does not
allocate temporary arrays.
NOTE 2: this algorithm uses BOTH previous points and current one, i.e.
new value of X[i] depends on BOTH previous point and X[i] itself.
NOTE 3: technical analytis users quite often work with EMA coefficient
expressed in DAYS instead of fractions. If you want to calculate
EMA(N), where N is a number of days, you can use alpha=2/(N+1).
-- ALGLIB --
Copyright 25.10.2011 by Bochkanov Sergey
*************************************************************************/
void filterema(real_1d_array &x, const ae_int_t n, const double alpha);
void filterema(real_1d_array &x, const double alpha);
/*************************************************************************
Filters: linear regression moving averages.
This filter replaces array by results of LRMA(K) filter.
LRMA(K) is defined as filter which, for each data point, builds linear
regression model using K prevous points (point itself is included in
these K points) and calculates value of this linear model at the point in
question.
INPUT PARAMETERS:
X - array[N], array to process. It can be larger than N,
in this case only first N points are processed.
N - points count, N>=0
K - K>=1 (K can be larger than N , such cases will be
correctly handled). Window width. K=1 corresponds to
identity transformation (nothing changes).
OUTPUT PARAMETERS:
X - array, whose first N elements were processed with SMA(K
)
NOTE 1: this function uses efficient in-place algorithm which does not
allocate temporary arrays.
NOTE 2: this algorithm makes only one pass through array and uses running
sum to speed-up calculation of the averages. Additional measures
are taken to ensure that running sum on a long sequence of zero
elements will be correctly reset to zero even in the presence of
round-off error.
NOTE 3: this is unsymmetric version of the algorithm, which does NOT
averages points after the current one. Only X[i], X[i-1], ... are
used when calculating new value of X[i]. We should also note that
this algorithm uses BOTH previous points and current one, i.e.
new value of X[i] depends on BOTH previous point and X[i] itself.
-- ALGLIB --
Copyright 25.10.2011 by Bochkanov Sergey
*************************************************************************/
void filterlrma(real_1d_array &x, const ae_int_t n, const ae_int_t k);
void filterlrma(real_1d_array &x, const ae_int_t k);
/*************************************************************************
k-means++ clusterization
INPUT PARAMETERS:
XY - dataset, array [0..NPoints-1,0..NVars-1].
NPoints - dataset size, NPoints>=K
NVars - number of variables, NVars>=1
K - desired number of clusters, K>=1
Restarts - number of restarts, Restarts>=1
OUTPUT PARAMETERS:
Info - return code:
* -3, if task is degenerate (number of distinct points
is
less than K)
* -1, if incorrect NPoints/NFeatures/K/Restarts was pas
sed
* 1, if subroutine finished successfully
C - array[0..NVars-1,0..K-1].matrix whose columns store
cluster's centers
XYC - array[NPoints], which contains cluster indexes
-- ALGLIB --
Copyright 21.03.2009 by Bochkanov Sergey
*************************************************************************/
void kmeansgenerate(const real_2d_array &xy, const ae_int_t npoints, const
ae_int_t nvars, const ae_int_t k, const ae_int_t restarts, ae_int_t &info,
real_2d_array &c, integer_1d_array &xyc);
/*************************************************************************
Multiclass Fisher LDA
Subroutine finds coefficients of linear combination which optimally separat
es
training set on classes.
INPUT PARAMETERS:
XY - training set, array[0..NPoints-1,0..NVars].
First NVars columns store values of independent
variables, next column stores number of class (from 0
to NClasses-1) which dataset element belongs to. Fracti
onal
values are rounded to nearest integer.
NPoints - training set size, NPoints>=0
NVars - number of independent variables, NVars>=1
NClasses - number of classes, NClasses>=2
OUTPUT PARAMETERS:
Info - return code:
* -4, if internal EVD subroutine hasn't converged
* -2, if there is a point with class number
outside of [0..NClasses-1].
* -1, if incorrect parameters was passed (NPoints<0,
NVars<1, NClasses<2)
* 1, if task has been solved
* 2, if there was a multicollinearity in training set,
but task has been solved.
W - linear combination coefficients, array[0..NVars-1]
-- ALGLIB --
Copyright 31.05.2008 by Bochkanov Sergey
*************************************************************************/
void fisherlda(const real_2d_array &xy, const ae_int_t npoints, const ae_in
t_t nvars, const ae_int_t nclasses, ae_int_t &info, real_1d_array &w);
/*************************************************************************
N-dimensional multiclass Fisher LDA
Subroutine finds coefficients of linear combinations which optimally separa
tes
training set on classes. It returns N-dimensional basis whose vector are so
rted
by quality of training set separation (in descending order).
INPUT PARAMETERS:
XY - training set, array[0..NPoints-1,0..NVars].
First NVars columns store values of independent
variables, next column stores number of class (from 0
to NClasses-1) which dataset element belongs to. Fracti
onal
values are rounded to nearest integer.
NPoints - training set size, NPoints>=0
NVars - number of independent variables, NVars>=1
NClasses - number of classes, NClasses>=2
OUTPUT PARAMETERS:
Info - return code:
* -4, if internal EVD subroutine hasn't converged
* -2, if there is a point with class number
outside of [0..NClasses-1].
* -1, if incorrect parameters was passed (NPoints<0,
NVars<1, NClasses<2)
* 1, if task has been solved
* 2, if there was a multicollinearity in training set,
but task has been solved.
W - basis, array[0..NVars-1,0..NVars-1]
columns of matrix stores basis vectors, sorted by
quality of training set separation (in descending order
)
-- ALGLIB --
Copyright 31.05.2008 by Bochkanov Sergey
*************************************************************************/
void fisherldan(const real_2d_array &xy, const ae_int_t npoints, const ae_i
nt_t nvars, const ae_int_t nclasses, ae_int_t &info, real_2d_array &w);
/*************************************************************************
This function serializes data structure to string. This function serializes data structure to string.
Important properties of s_out: Important properties of s_out:
* it contains alphanumeric characters, dots, underscores, minus signs * it contains alphanumeric characters, dots, underscores, minus signs
* these symbols are grouped into words, which are separated by spaces * these symbols are grouped into words, which are separated by spaces
and Windows-style (CR+LF) newlines and Windows-style (CR+LF) newlines
* although serializer uses spaces and CR+LF as separators, you can * although serializer uses spaces and CR+LF as separators, you can
replace any separator character by arbitrary combination of spaces, replace any separator character by arbitrary combination of spaces,
tabs, Windows or Unix newlines. It allows flexible reformatting of tabs, Windows or Unix newlines. It allows flexible reformatting of
the string in case you want to include it into text or XML file. the string in case you want to include it into text or XML file.
skipping to change at line 1332 skipping to change at line 1435
/************************************************************************* /*************************************************************************
Returns information about initialized network: number of inputs, outputs, Returns information about initialized network: number of inputs, outputs,
weights. weights.
-- ALGLIB -- -- ALGLIB --
Copyright 04.11.2007 by Bochkanov Sergey Copyright 04.11.2007 by Bochkanov Sergey
*************************************************************************/ *************************************************************************/
void mlpproperties(const multilayerperceptron &network, ae_int_t &nin, ae_i nt_t &nout, ae_int_t &wcount); void mlpproperties(const multilayerperceptron &network, ae_int_t &nin, ae_i nt_t &nout, ae_int_t &wcount);
/************************************************************************* /*************************************************************************
Returns number of inputs.
-- ALGLIB --
Copyright 19.10.2011 by Bochkanov Sergey
*************************************************************************/
ae_int_t mlpgetinputscount(const multilayerperceptron &network);
/*************************************************************************
Returns number of outputs.
-- ALGLIB --
Copyright 19.10.2011 by Bochkanov Sergey
*************************************************************************/
ae_int_t mlpgetoutputscount(const multilayerperceptron &network);
/*************************************************************************
Returns number of weights.
-- ALGLIB --
Copyright 19.10.2011 by Bochkanov Sergey
*************************************************************************/
ae_int_t mlpgetweightscount(const multilayerperceptron &network);
/*************************************************************************
Tells whether network is SOFTMAX-normalized (i.e. classifier) or not. Tells whether network is SOFTMAX-normalized (i.e. classifier) or not.
-- ALGLIB -- -- ALGLIB --
Copyright 04.11.2007 by Bochkanov Sergey Copyright 04.11.2007 by Bochkanov Sergey
*************************************************************************/ *************************************************************************/
bool mlpissoftmax(const multilayerperceptron &network); bool mlpissoftmax(const multilayerperceptron &network);
/************************************************************************* /*************************************************************************
This function returns total number of layers (including input, hidden and This function returns total number of layers (including input, hidden and
output layers). output layers).
skipping to change at line 2700 skipping to change at line 2827
*************************************************************************/ *************************************************************************/
void mlptrainlbfgs(const multilayerperceptron &network, const real_2d_array &xy, const ae_int_t npoints, const double decay, const ae_int_t restarts, const double wstep, const ae_int_t maxits, ae_int_t &info, mlpreport &rep); void mlptrainlbfgs(const multilayerperceptron &network, const real_2d_array &xy, const ae_int_t npoints, const double decay, const ae_int_t restarts, const double wstep, const ae_int_t maxits, ae_int_t &info, mlpreport &rep);
/************************************************************************* /*************************************************************************
Neural network training using early stopping (base algorithm - L-BFGS with Neural network training using early stopping (base algorithm - L-BFGS with
regularization). regularization).
INPUT PARAMETERS: INPUT PARAMETERS:
Network - neural network with initialized geometry Network - neural network with initialized geometry
TrnXY - training set TrnXY - training set
TrnSize - training set size TrnSize - training set size, TrnSize>0
ValXY - validation set ValXY - validation set
ValSize - validation set size ValSize - validation set size, ValSize>0
Decay - weight decay constant, >=0.001 Decay - weight decay constant, >=0.001
Decay term 'Decay*||Weights||^2' is added to error Decay term 'Decay*||Weights||^2' is added to error
function. function.
If you don't know what Decay to choose, use 0.001. If you don't know what Decay to choose, use 0.001.
Restarts - number of restarts from random position, >0. Restarts - number of restarts, either:
If you don't know what Restarts to choose, use 2. * strictly positive number - algorithm make specified
number of restarts from random position.
* -1, in which case algorithm makes exactly one run
from the initial state of the network (no randomizati
on).
If you don't know what Restarts to choose, choose one
one the following:
* -1 (deterministic start)
* +1 (one random restart)
* +5 (moderate amount of random restarts)
OUTPUT PARAMETERS: OUTPUT PARAMETERS:
Network - trained neural network. Network - trained neural network.
Info - return code: Info - return code:
* -2, if there is a point with class number * -2, if there is a point with class number
outside of [0..NOut-1]. outside of [0..NOut-1].
* -1, if wrong parameters specified * -1, if wrong parameters specified
(NPoints<0, Restarts<1, ...). (NPoints<0, Restarts<1, ...).
* 2, task has been solved, stopping criterion met - * 2, task has been solved, stopping criterion met -
sufficiently small step size. Not expected (we sufficiently small step size. Not expected (we
skipping to change at line 2798 skipping to change at line 2933
Info - return code, same as in MLPTrainLBFGS Info - return code, same as in MLPTrainLBFGS
Rep - report, same as in MLPTrainLM/MLPTrainLBFGS Rep - report, same as in MLPTrainLM/MLPTrainLBFGS
CVRep - generalization error estimates CVRep - generalization error estimates
-- ALGLIB -- -- ALGLIB --
Copyright 09.12.2007 by Bochkanov Sergey Copyright 09.12.2007 by Bochkanov Sergey
*************************************************************************/ *************************************************************************/
void mlpkfoldcvlm(const multilayerperceptron &network, const real_2d_array &xy, const ae_int_t npoints, const double decay, const ae_int_t restarts, c onst ae_int_t foldscount, ae_int_t &info, mlpreport &rep, mlpcvreport &cvre p); void mlpkfoldcvlm(const multilayerperceptron &network, const real_2d_array &xy, const ae_int_t npoints, const double decay, const ae_int_t restarts, c onst ae_int_t foldscount, ae_int_t &info, mlpreport &rep, mlpcvreport &cvre p);
/************************************************************************* /*************************************************************************
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.
*************************************************************************/
void mlpeserialize(mlpensemble &obj, std::string &s_out);
/*************************************************************************
This function unserializes data structure from string.
*************************************************************************/
void mlpeunserialize(std::string &s_in, mlpensemble &obj);
/*************************************************************************
Like MLPCreate0, but for ensembles. Like MLPCreate0, but for ensembles.
-- ALGLIB -- -- ALGLIB --
Copyright 18.02.2009 by Bochkanov Sergey Copyright 18.02.2009 by Bochkanov Sergey
*************************************************************************/ *************************************************************************/
void mlpecreate0(const ae_int_t nin, const ae_int_t nout, const ae_int_t en semblesize, mlpensemble &ensemble); void mlpecreate0(const ae_int_t nin, const ae_int_t nout, const ae_int_t en semblesize, mlpensemble &ensemble);
/************************************************************************* /*************************************************************************
Like MLPCreate1, but for ensembles. Like MLPCreate1, but for ensembles.
skipping to change at line 3334 skipping to change at line 3496
ae_state *_state); ae_state *_state);
ae_bool _decisionforest_init(decisionforest* p, ae_state *_state, ae_bool m ake_automatic); ae_bool _decisionforest_init(decisionforest* p, ae_state *_state, ae_bool m ake_automatic);
ae_bool _decisionforest_init_copy(decisionforest* dst, decisionforest* src, ae_state *_state, ae_bool make_automatic); ae_bool _decisionforest_init_copy(decisionforest* dst, decisionforest* src, ae_state *_state, ae_bool make_automatic);
void _decisionforest_clear(decisionforest* p); void _decisionforest_clear(decisionforest* p);
ae_bool _dfreport_init(dfreport* p, ae_state *_state, ae_bool make_automati c); ae_bool _dfreport_init(dfreport* p, ae_state *_state, ae_bool make_automati c);
ae_bool _dfreport_init_copy(dfreport* dst, dfreport* src, ae_state *_state, ae_bool make_automatic); ae_bool _dfreport_init_copy(dfreport* dst, dfreport* src, ae_state *_state, ae_bool make_automatic);
void _dfreport_clear(dfreport* p); void _dfreport_clear(dfreport* p);
ae_bool _dfinternalbuffers_init(dfinternalbuffers* p, ae_state *_state, ae_ bool make_automatic); ae_bool _dfinternalbuffers_init(dfinternalbuffers* p, ae_state *_state, ae_ bool make_automatic);
ae_bool _dfinternalbuffers_init_copy(dfinternalbuffers* dst, dfinternalbuff ers* src, ae_state *_state, ae_bool make_automatic); ae_bool _dfinternalbuffers_init_copy(dfinternalbuffers* dst, dfinternalbuff ers* src, ae_state *_state, ae_bool make_automatic);
void _dfinternalbuffers_clear(dfinternalbuffers* p); void _dfinternalbuffers_clear(dfinternalbuffers* p);
void kmeansgenerate(/* Real */ ae_matrix* xy,
ae_int_t npoints,
ae_int_t nvars,
ae_int_t k,
ae_int_t restarts,
ae_int_t* info,
/* Real */ ae_matrix* c,
/* Integer */ ae_vector* xyc,
ae_state *_state);
void fisherlda(/* Real */ ae_matrix* xy,
ae_int_t npoints,
ae_int_t nvars,
ae_int_t nclasses,
ae_int_t* info,
/* Real */ ae_vector* w,
ae_state *_state);
void fisherldan(/* Real */ ae_matrix* xy,
ae_int_t npoints,
ae_int_t nvars,
ae_int_t nclasses,
ae_int_t* info,
/* Real */ ae_matrix* w,
ae_state *_state);
void lrbuild(/* Real */ ae_matrix* xy, void lrbuild(/* Real */ ae_matrix* xy,
ae_int_t npoints, ae_int_t npoints,
ae_int_t nvars, ae_int_t nvars,
ae_int_t* info, ae_int_t* info,
linearmodel* lm, linearmodel* lm,
lrreport* ar, lrreport* ar,
ae_state *_state); ae_state *_state);
void lrbuilds(/* Real */ ae_matrix* xy, void lrbuilds(/* Real */ ae_matrix* xy,
/* Real */ ae_vector* s, /* Real */ ae_vector* s,
ae_int_t npoints, ae_int_t npoints,
skipping to change at line 3435 skipping to change at line 3574
ae_int_t* info, ae_int_t* info,
double* a, double* a,
double* b, double* b,
ae_state *_state); ae_state *_state);
ae_bool _linearmodel_init(linearmodel* p, ae_state *_state, ae_bool make_au tomatic); ae_bool _linearmodel_init(linearmodel* p, ae_state *_state, ae_bool make_au tomatic);
ae_bool _linearmodel_init_copy(linearmodel* dst, linearmodel* src, ae_state *_state, ae_bool make_automatic); ae_bool _linearmodel_init_copy(linearmodel* dst, linearmodel* src, ae_state *_state, ae_bool make_automatic);
void _linearmodel_clear(linearmodel* p); void _linearmodel_clear(linearmodel* p);
ae_bool _lrreport_init(lrreport* p, ae_state *_state, ae_bool make_automati c); ae_bool _lrreport_init(lrreport* p, ae_state *_state, ae_bool make_automati c);
ae_bool _lrreport_init_copy(lrreport* dst, lrreport* src, ae_state *_state, ae_bool make_automatic); ae_bool _lrreport_init_copy(lrreport* dst, lrreport* src, ae_state *_state, ae_bool make_automatic);
void _lrreport_clear(lrreport* p); void _lrreport_clear(lrreport* p);
void filtersma(/* Real */ ae_vector* x,
ae_int_t n,
ae_int_t k,
ae_state *_state);
void filterema(/* Real */ ae_vector* x,
ae_int_t n,
double alpha,
ae_state *_state);
void filterlrma(/* Real */ ae_vector* x,
ae_int_t n,
ae_int_t k,
ae_state *_state);
void kmeansgenerate(/* Real */ ae_matrix* xy,
ae_int_t npoints,
ae_int_t nvars,
ae_int_t k,
ae_int_t restarts,
ae_int_t* info,
/* Real */ ae_matrix* c,
/* Integer */ ae_vector* xyc,
ae_state *_state);
void fisherlda(/* Real */ ae_matrix* xy,
ae_int_t npoints,
ae_int_t nvars,
ae_int_t nclasses,
ae_int_t* info,
/* Real */ ae_vector* w,
ae_state *_state);
void fisherldan(/* Real */ ae_matrix* xy,
ae_int_t npoints,
ae_int_t nvars,
ae_int_t nclasses,
ae_int_t* info,
/* Real */ ae_matrix* w,
ae_state *_state);
void mlpcreate0(ae_int_t nin, void mlpcreate0(ae_int_t nin,
ae_int_t nout, ae_int_t nout,
multilayerperceptron* network, multilayerperceptron* network,
ae_state *_state); ae_state *_state);
void mlpcreate1(ae_int_t nin, void mlpcreate1(ae_int_t nin,
ae_int_t nhid, ae_int_t nhid,
ae_int_t nout, ae_int_t nout,
multilayerperceptron* network, multilayerperceptron* network,
ae_state *_state); ae_state *_state);
void mlpcreate2(ae_int_t nin, void mlpcreate2(ae_int_t nin,
skipping to change at line 3528 skipping to change at line 3702
void mlprandomizefull(multilayerperceptron* network, ae_state *_state); void mlprandomizefull(multilayerperceptron* network, ae_state *_state);
void mlpinitpreprocessor(multilayerperceptron* network, void mlpinitpreprocessor(multilayerperceptron* network,
/* Real */ ae_matrix* xy, /* Real */ ae_matrix* xy,
ae_int_t ssize, ae_int_t ssize,
ae_state *_state); ae_state *_state);
void mlpproperties(multilayerperceptron* network, void mlpproperties(multilayerperceptron* network,
ae_int_t* nin, ae_int_t* nin,
ae_int_t* nout, ae_int_t* nout,
ae_int_t* wcount, ae_int_t* wcount,
ae_state *_state); ae_state *_state);
ae_int_t mlpgetinputscount(multilayerperceptron* network,
ae_state *_state);
ae_int_t mlpgetoutputscount(multilayerperceptron* network,
ae_state *_state);
ae_int_t mlpgetweightscount(multilayerperceptron* network,
ae_state *_state);
ae_bool mlpissoftmax(multilayerperceptron* network, ae_state *_state); ae_bool mlpissoftmax(multilayerperceptron* network, ae_state *_state);
ae_int_t mlpgetlayerscount(multilayerperceptron* network, ae_int_t mlpgetlayerscount(multilayerperceptron* network,
ae_state *_state); ae_state *_state);
ae_int_t mlpgetlayersize(multilayerperceptron* network, ae_int_t mlpgetlayersize(multilayerperceptron* network,
ae_int_t k, ae_int_t k,
ae_state *_state); ae_state *_state);
void mlpgetinputscaling(multilayerperceptron* network, void mlpgetinputscaling(multilayerperceptron* network,
ae_int_t i, ae_int_t i,
double* mean, double* mean,
double* sigma, double* sigma,
skipping to change at line 3950 skipping to change at line 4130
ae_int_t ensemblesize, ae_int_t ensemblesize,
mlpensemble* ensemble, mlpensemble* ensemble,
ae_state *_state); ae_state *_state);
void mlpecreatefromnetwork(multilayerperceptron* network, void mlpecreatefromnetwork(multilayerperceptron* network,
ae_int_t ensemblesize, ae_int_t ensemblesize,
mlpensemble* ensemble, mlpensemble* ensemble,
ae_state *_state); ae_state *_state);
void mlpecopy(mlpensemble* ensemble1, void mlpecopy(mlpensemble* ensemble1,
mlpensemble* ensemble2, mlpensemble* ensemble2,
ae_state *_state); ae_state *_state);
void mlpeserialize(mlpensemble* ensemble,
/* Real */ ae_vector* ra,
ae_int_t* rlen,
ae_state *_state);
void mlpeunserialize(/* Real */ ae_vector* ra,
mlpensemble* ensemble,
ae_state *_state);
void mlperandomize(mlpensemble* ensemble, ae_state *_state); void mlperandomize(mlpensemble* ensemble, ae_state *_state);
void mlpeproperties(mlpensemble* ensemble, void mlpeproperties(mlpensemble* ensemble,
ae_int_t* nin, ae_int_t* nin,
ae_int_t* nout, ae_int_t* nout,
ae_state *_state); ae_state *_state);
ae_bool mlpeissoftmax(mlpensemble* ensemble, ae_state *_state); ae_bool mlpeissoftmax(mlpensemble* ensemble, ae_state *_state);
void mlpeprocess(mlpensemble* ensemble, void mlpeprocess(mlpensemble* ensemble,
/* Real */ ae_vector* x, /* Real */ ae_vector* x,
/* Real */ ae_vector* y, /* Real */ ae_vector* y,
ae_state *_state); ae_state *_state);
skipping to change at line 4019 skipping to change at line 4192
mlpcvreport* ooberrors, mlpcvreport* ooberrors,
ae_state *_state); ae_state *_state);
void mlpetraines(mlpensemble* ensemble, void mlpetraines(mlpensemble* ensemble,
/* Real */ ae_matrix* xy, /* Real */ ae_matrix* xy,
ae_int_t npoints, ae_int_t npoints,
double decay, double decay,
ae_int_t restarts, ae_int_t restarts,
ae_int_t* info, ae_int_t* info,
mlpreport* rep, mlpreport* rep,
ae_state *_state); ae_state *_state);
void mlpealloc(ae_serializer* s, mlpensemble* ensemble, ae_state *_state);
void mlpeserialize(ae_serializer* s,
mlpensemble* ensemble,
ae_state *_state);
void mlpeunserialize(ae_serializer* s,
mlpensemble* ensemble,
ae_state *_state);
ae_bool _mlpensemble_init(mlpensemble* p, ae_state *_state, ae_bool make_au tomatic); ae_bool _mlpensemble_init(mlpensemble* p, ae_state *_state, ae_bool make_au tomatic);
ae_bool _mlpensemble_init_copy(mlpensemble* dst, mlpensemble* src, ae_state *_state, ae_bool make_automatic); ae_bool _mlpensemble_init_copy(mlpensemble* dst, mlpensemble* src, ae_state *_state, ae_bool make_automatic);
void _mlpensemble_clear(mlpensemble* p); void _mlpensemble_clear(mlpensemble* p);
void pcabuildbasis(/* Real */ ae_matrix* x, void pcabuildbasis(/* Real */ ae_matrix* x,
ae_int_t npoints, ae_int_t npoints,
ae_int_t nvars, ae_int_t nvars,
ae_int_t* info, ae_int_t* info,
/* Real */ ae_vector* s2, /* Real */ ae_vector* s2,
/* Real */ ae_matrix* v, /* Real */ ae_matrix* v,
ae_state *_state); ae_state *_state);
 End of changes. 17 change blocks. 
154 lines changed or deleted 337 lines changed or added


 interpolation.h   interpolation.h 
skipping to change at line 68 skipping to change at line 68
double sy; double sy;
ae_vector x; ae_vector x;
ae_vector y; ae_vector y;
ae_vector w; ae_vector w;
} barycentricinterpolant; } barycentricinterpolant;
typedef struct typedef struct
{ {
ae_bool periodic; ae_bool periodic;
ae_int_t n; ae_int_t n;
ae_int_t k; ae_int_t k;
ae_int_t continuity;
ae_vector x; ae_vector x;
ae_vector c; ae_vector c;
} spline1dinterpolant; } spline1dinterpolant;
typedef struct typedef struct
{ {
double taskrcond; double taskrcond;
double rmserror; double rmserror;
double avgerror; double avgerror;
double avgrelerror; double avgrelerror;
double maxerror; double maxerror;
skipping to change at line 168 skipping to change at line 169
{ {
ae_int_t n; ae_int_t n;
ae_bool periodic; ae_bool periodic;
ae_vector p; ae_vector p;
spline1dinterpolant x; spline1dinterpolant x;
spline1dinterpolant y; spline1dinterpolant y;
spline1dinterpolant z; spline1dinterpolant z;
} pspline3interpolant; } pspline3interpolant;
typedef struct typedef struct
{ {
ae_int_t ny;
ae_int_t nx;
ae_int_t nc;
ae_int_t nl;
kdtree tree;
ae_matrix xc;
ae_matrix wr;
double rmax;
ae_matrix v;
ae_int_t gridtype;
ae_bool fixrad;
double lambdav;
double radvalue;
double radzvalue;
ae_int_t nlayers;
ae_int_t aterm;
ae_int_t algorithmtype;
double epsort;
double epserr;
ae_int_t maxits;
double h;
ae_int_t n;
ae_matrix x;
ae_matrix y;
ae_vector calcbufxcx;
ae_matrix calcbufx;
ae_vector calcbuftags;
} rbfmodel;
typedef struct
{
ae_int_t arows;
ae_int_t acols;
ae_int_t annz;
ae_int_t iterationscount;
ae_int_t nmv;
ae_int_t terminationtype;
} rbfreport;
typedef struct
{
ae_int_t k; ae_int_t k;
ae_vector c; ae_vector c;
} spline2dinterpolant; } spline2dinterpolant;
} }
///////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////
// //
// THIS SECTION CONTAINS C++ INTERFACE // THIS SECTION CONTAINS C++ INTERFACE
// //
skipping to change at line 233 skipping to change at line 273
{ {
public: public:
barycentricinterpolant(); barycentricinterpolant();
barycentricinterpolant(const barycentricinterpolant &rhs); barycentricinterpolant(const barycentricinterpolant &rhs);
barycentricinterpolant& operator=(const barycentricinterpolant &rhs); barycentricinterpolant& operator=(const barycentricinterpolant &rhs);
virtual ~barycentricinterpolant(); virtual ~barycentricinterpolant();
}; };
/************************************************************************* /*************************************************************************
1-dimensional spline inteprolant 1-dimensional spline interpolant
*************************************************************************/ *************************************************************************/
class _spline1dinterpolant_owner class _spline1dinterpolant_owner
{ {
public: public:
_spline1dinterpolant_owner(); _spline1dinterpolant_owner();
_spline1dinterpolant_owner(const _spline1dinterpolant_owner &rhs); _spline1dinterpolant_owner(const _spline1dinterpolant_owner &rhs);
_spline1dinterpolant_owner& operator=(const _spline1dinterpolant_owner &rhs); _spline1dinterpolant_owner& operator=(const _spline1dinterpolant_owner &rhs);
virtual ~_spline1dinterpolant_owner(); virtual ~_spline1dinterpolant_owner();
alglib_impl::spline1dinterpolant* c_ptr(); alglib_impl::spline1dinterpolant* c_ptr();
alglib_impl::spline1dinterpolant* c_ptr() const; alglib_impl::spline1dinterpolant* c_ptr() const;
skipping to change at line 501 skipping to change at line 541
{ {
public: public:
pspline3interpolant(); pspline3interpolant();
pspline3interpolant(const pspline3interpolant &rhs); pspline3interpolant(const pspline3interpolant &rhs);
pspline3interpolant& operator=(const pspline3interpolant &rhs); pspline3interpolant& operator=(const pspline3interpolant &rhs);
virtual ~pspline3interpolant(); virtual ~pspline3interpolant();
}; };
/************************************************************************* /*************************************************************************
RBF model.
Never try to directly work with fields of this object - always use ALGLIB
functions to use this object.
*************************************************************************/
class _rbfmodel_owner
{
public:
_rbfmodel_owner();
_rbfmodel_owner(const _rbfmodel_owner &rhs);
_rbfmodel_owner& operator=(const _rbfmodel_owner &rhs);
virtual ~_rbfmodel_owner();
alglib_impl::rbfmodel* c_ptr();
alglib_impl::rbfmodel* c_ptr() const;
protected:
alglib_impl::rbfmodel *p_struct;
};
class rbfmodel : public _rbfmodel_owner
{
public:
rbfmodel();
rbfmodel(const rbfmodel &rhs);
rbfmodel& operator=(const rbfmodel &rhs);
virtual ~rbfmodel();
};
/*************************************************************************
RBF solution report:
* TerminationType - termination type, positive values - success,
non-positive - failure.
*************************************************************************/
class _rbfreport_owner
{
public:
_rbfreport_owner();
_rbfreport_owner(const _rbfreport_owner &rhs);
_rbfreport_owner& operator=(const _rbfreport_owner &rhs);
virtual ~_rbfreport_owner();
alglib_impl::rbfreport* c_ptr();
alglib_impl::rbfreport* c_ptr() const;
protected:
alglib_impl::rbfreport *p_struct;
};
class rbfreport : public _rbfreport_owner
{
public:
rbfreport();
rbfreport(const rbfreport &rhs);
rbfreport& operator=(const rbfreport &rhs);
virtual ~rbfreport();
ae_int_t &arows;
ae_int_t &acols;
ae_int_t &annz;
ae_int_t &iterationscount;
ae_int_t &nmv;
ae_int_t &terminationtype;
};
/*************************************************************************
2-dimensional spline inteprolant 2-dimensional spline inteprolant
*************************************************************************/ *************************************************************************/
class _spline2dinterpolant_owner class _spline2dinterpolant_owner
{ {
public: public:
_spline2dinterpolant_owner(); _spline2dinterpolant_owner();
_spline2dinterpolant_owner(const _spline2dinterpolant_owner &rhs); _spline2dinterpolant_owner(const _spline2dinterpolant_owner &rhs);
_spline2dinterpolant_owner& operator=(const _spline2dinterpolant_owner &rhs); _spline2dinterpolant_owner& operator=(const _spline2dinterpolant_owner &rhs);
virtual ~_spline2dinterpolant_owner(); virtual ~_spline2dinterpolant_owner();
alglib_impl::spline2dinterpolant* c_ptr(); alglib_impl::spline2dinterpolant* c_ptr();
skipping to change at line 3358 skipping to change at line 3459
RESULT: RESULT:
length of arc starting at T=A and ending at T=B. length of arc starting at T=A and ending at T=B.
-- ALGLIB PROJECT -- -- ALGLIB PROJECT --
Copyright 30.05.2010 by Bochkanov Sergey Copyright 30.05.2010 by Bochkanov Sergey
*************************************************************************/ *************************************************************************/
double pspline3arclength(const pspline3interpolant &p, const double a, cons t double b); double pspline3arclength(const pspline3interpolant &p, const double a, cons t double b);
/************************************************************************* /*************************************************************************
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.
*************************************************************************/
void rbfserialize(rbfmodel &obj, std::string &s_out);
/*************************************************************************
This function unserializes data structure from string.
*************************************************************************/
void rbfunserialize(std::string &s_in, rbfmodel &obj);
/*************************************************************************
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
*************************************************************************/
void rbfcreate(const ae_int_t nx, const ae_int_t ny, rbfmodel &s);
/*************************************************************************
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
*************************************************************************/
void rbfsetpoints(const rbfmodel &s, const real_2d_array &xy, const ae_int_
t n);
void rbfsetpoints(const rbfmodel &s, const real_2d_array &xy);
/*************************************************************************
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
*************************************************************************/
void rbfsetalgoqnn(const rbfmodel &s, const double q, const double z);
void rbfsetalgoqnn(const rbfmodel &s);
/*************************************************************************
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. Using too many layers - you will get too large and too slow model. 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
*************************************************************************/
void rbfsetalgomultilayer(const rbfmodel &s, const double rbase, const ae_i
nt_t nlayers, const double lambdav);
void rbfsetalgomultilayer(const rbfmodel &s, const double rbase, const ae_i
nt_t nlayers);
/*************************************************************************
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
*************************************************************************/
void rbfsetlinterm(const rbfmodel &s);
/*************************************************************************
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
*************************************************************************/
void rbfsetconstterm(const rbfmodel &s);
/*************************************************************************
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
*************************************************************************/
void rbfsetzeroterm(const rbfmodel &s);
/*************************************************************************
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
*************************************************************************/
void rbfbuildmodel(const rbfmodel &s, rbfreport &rep);
/*************************************************************************
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
*************************************************************************/
double rbfcalc2(const rbfmodel &s, const double x0, const double x1);
/*************************************************************************
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
*************************************************************************/
double rbfcalc3(const rbfmodel &s, const double x0, const double x1, const
double x2);
/*************************************************************************
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
*************************************************************************/
void rbfcalc(const rbfmodel &s, const real_1d_array &x, real_1d_array &y);
/*************************************************************************
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
*************************************************************************/
void rbfcalcbuf(const rbfmodel &s, const real_1d_array &x, real_1d_array &y
);
/*************************************************************************
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
*************************************************************************/
void rbfgridcalc2(const rbfmodel &s, const real_1d_array &x0, const ae_int_
t n0, const real_1d_array &x1, const ae_int_t n1, real_2d_array &y);
/*************************************************************************
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
*************************************************************************/
void rbfunpack(const rbfmodel &s, ae_int_t &nx, ae_int_t &ny, real_2d_array
&xwr, ae_int_t &nc, real_2d_array &v);
/*************************************************************************
This subroutine builds bilinear spline coefficients table. This subroutine builds bilinear spline coefficients table.
Input parameters: Input parameters:
X - spline abscissas, array[0..N-1] X - spline abscissas, array[0..N-1]
Y - spline ordinates, array[0..M-1] Y - spline ordinates, array[0..M-1]
F - function values, array[0..M-1,0..N-1] F - function values, array[0..M-1,0..N-1]
M,N - grid size, M>=2, N>=2 M,N - grid size, M>=2, N>=2
Output parameters: Output parameters:
C - spline interpolant C - spline interpolant
skipping to change at line 3803 skipping to change at line 4404
ae_int_t n, ae_int_t n,
/* Real */ ae_vector* x2, /* Real */ ae_vector* x2,
ae_int_t n2, ae_int_t n2,
/* Real */ ae_vector* y, /* Real */ ae_vector* y,
ae_bool needy, ae_bool needy,
/* Real */ ae_vector* d1, /* Real */ ae_vector* d1,
ae_bool needd1, ae_bool needd1,
/* Real */ ae_vector* d2, /* Real */ ae_vector* d2,
ae_bool needd2, ae_bool needd2,
ae_state *_state); ae_state *_state);
void spline1drootsandextrema(spline1dinterpolant* c,
/* Real */ ae_vector* r,
ae_int_t* nr,
ae_bool* dr,
/* Real */ ae_vector* e,
/* Integer */ ae_vector* et,
ae_int_t* ne,
ae_bool* de,
ae_state *_state);
void heapsortdpoints(/* Real */ ae_vector* x, void heapsortdpoints(/* Real */ ae_vector* x,
/* Real */ ae_vector* y, /* Real */ ae_vector* y,
/* Real */ ae_vector* d, /* Real */ ae_vector* d,
ae_int_t n, ae_int_t n,
ae_state *_state); ae_state *_state);
void solvepolinom2(double p0,
double m0,
double p1,
double m1,
double* x0,
double* x1,
ae_int_t* nr,
ae_state *_state);
void solvecubicpolinom(double pa,
double ma,
double pb,
double mb,
double a,
double b,
double* x0,
double* x1,
double* x2,
double* ex0,
double* ex1,
ae_int_t* nr,
ae_int_t* ne,
/* Real */ ae_vector* tempdata,
ae_state *_state);
ae_int_t bisectmethod(double pa,
double ma,
double pb,
double mb,
double a,
double b,
double* x,
ae_state *_state);
ae_bool _spline1dinterpolant_init(spline1dinterpolant* p, ae_state *_state, ae_bool make_automatic); ae_bool _spline1dinterpolant_init(spline1dinterpolant* p, ae_state *_state, ae_bool make_automatic);
ae_bool _spline1dinterpolant_init_copy(spline1dinterpolant* dst, spline1din terpolant* src, ae_state *_state, ae_bool make_automatic); ae_bool _spline1dinterpolant_init_copy(spline1dinterpolant* dst, spline1din terpolant* src, ae_state *_state, ae_bool make_automatic);
void _spline1dinterpolant_clear(spline1dinterpolant* p); void _spline1dinterpolant_clear(spline1dinterpolant* p);
void polynomialfit(/* Real */ ae_vector* x, void polynomialfit(/* Real */ ae_vector* x,
/* Real */ ae_vector* y, /* Real */ ae_vector* y,
ae_int_t n, ae_int_t n,
ae_int_t m, ae_int_t m,
ae_int_t* info, ae_int_t* info,
barycentricinterpolant* p, barycentricinterpolant* p,
polynomialfitreport* rep, polynomialfitreport* rep,
skipping to change at line 4162 skipping to change at line 4803
double pspline3arclength(pspline3interpolant* p, double pspline3arclength(pspline3interpolant* p,
double a, double a,
double b, double b,
ae_state *_state); ae_state *_state);
ae_bool _pspline2interpolant_init(pspline2interpolant* p, ae_state *_state, ae_bool make_automatic); ae_bool _pspline2interpolant_init(pspline2interpolant* p, ae_state *_state, ae_bool make_automatic);
ae_bool _pspline2interpolant_init_copy(pspline2interpolant* dst, pspline2in terpolant* src, ae_state *_state, ae_bool make_automatic); ae_bool _pspline2interpolant_init_copy(pspline2interpolant* dst, pspline2in terpolant* src, ae_state *_state, ae_bool make_automatic);
void _pspline2interpolant_clear(pspline2interpolant* p); void _pspline2interpolant_clear(pspline2interpolant* p);
ae_bool _pspline3interpolant_init(pspline3interpolant* p, ae_state *_state, ae_bool make_automatic); ae_bool _pspline3interpolant_init(pspline3interpolant* p, ae_state *_state, ae_bool make_automatic);
ae_bool _pspline3interpolant_init_copy(pspline3interpolant* dst, pspline3in terpolant* src, ae_state *_state, ae_bool make_automatic); ae_bool _pspline3interpolant_init_copy(pspline3interpolant* dst, pspline3in terpolant* src, ae_state *_state, ae_bool make_automatic);
void _pspline3interpolant_clear(pspline3interpolant* p); void _pspline3interpolant_clear(pspline3interpolant* p);
void rbfcreate(ae_int_t nx, ae_int_t ny, rbfmodel* s, ae_state *_state);
void rbfsetpoints(rbfmodel* s,
/* Real */ ae_matrix* xy,
ae_int_t n,
ae_state *_state);
void rbfsetalgoqnn(rbfmodel* s, double q, double z, ae_state *_state);
void rbfsetalgomultilayer(rbfmodel* s,
double rbase,
ae_int_t nlayers,
double lambdav,
ae_state *_state);
void rbfsetlinterm(rbfmodel* s, ae_state *_state);
void rbfsetconstterm(rbfmodel* s, ae_state *_state);
void rbfsetzeroterm(rbfmodel* s, ae_state *_state);
void rbfsetcond(rbfmodel* s,
double epsort,
double epserr,
ae_int_t maxits,
ae_state *_state);
void rbfbuildmodel(rbfmodel* s, rbfreport* rep, ae_state *_state);
double rbfcalc2(rbfmodel* s, double x0, double x1, ae_state *_state);
double rbfcalc3(rbfmodel* s,
double x0,
double x1,
double x2,
ae_state *_state);
void rbfcalc(rbfmodel* s,
/* Real */ ae_vector* x,
/* Real */ ae_vector* y,
ae_state *_state);
void rbfcalcbuf(rbfmodel* s,
/* Real */ ae_vector* x,
/* Real */ ae_vector* y,
ae_state *_state);
void rbfgridcalc2(rbfmodel* s,
/* Real */ ae_vector* x0,
ae_int_t n0,
/* Real */ ae_vector* x1,
ae_int_t n1,
/* Real */ ae_matrix* y,
ae_state *_state);
void rbfunpack(rbfmodel* s,
ae_int_t* nx,
ae_int_t* ny,
/* Real */ ae_matrix* xwr,
ae_int_t* nc,
/* Real */ ae_matrix* v,
ae_state *_state);
void rbfalloc(ae_serializer* s, rbfmodel* model, ae_state *_state);
void rbfserialize(ae_serializer* s, rbfmodel* model, ae_state *_state);
void rbfunserialize(ae_serializer* s, rbfmodel* model, ae_state *_state);
ae_bool _rbfmodel_init(rbfmodel* p, ae_state *_state, ae_bool make_automati
c);
ae_bool _rbfmodel_init_copy(rbfmodel* dst, rbfmodel* src, ae_state *_state,
ae_bool make_automatic);
void _rbfmodel_clear(rbfmodel* p);
ae_bool _rbfreport_init(rbfreport* p, ae_state *_state, ae_bool make_automa
tic);
ae_bool _rbfreport_init_copy(rbfreport* dst, rbfreport* src, ae_state *_sta
te, ae_bool make_automatic);
void _rbfreport_clear(rbfreport* p);
void spline2dbuildbilinear(/* Real */ ae_vector* x, void spline2dbuildbilinear(/* Real */ ae_vector* x,
/* Real */ ae_vector* y, /* Real */ ae_vector* y,
/* Real */ ae_matrix* f, /* Real */ ae_matrix* f,
ae_int_t m, ae_int_t m,
ae_int_t n, ae_int_t n,
spline2dinterpolant* c, spline2dinterpolant* c,
ae_state *_state); ae_state *_state);
void spline2dbuildbicubic(/* Real */ ae_vector* x, void spline2dbuildbicubic(/* Real */ ae_vector* x,
/* Real */ ae_vector* y, /* Real */ ae_vector* y,
/* Real */ ae_matrix* f, /* Real */ ae_matrix* f,
 End of changes. 8 change blocks. 
1 lines changed or deleted 710 lines changed or added


 linalg.h   linalg.h 
skipping to change at line 55 skipping to change at line 55
ae_vector rk; ae_vector rk;
ae_vector rk1; ae_vector rk1;
ae_vector xk; ae_vector xk;
ae_vector xk1; ae_vector xk1;
ae_vector pk; ae_vector pk;
ae_vector pk1; ae_vector pk1;
ae_vector b; ae_vector b;
rcommstate rstate; rcommstate rstate;
ae_vector tmp2; ae_vector tmp2;
} fblslincgstate; } fblslincgstate;
typedef struct
{
ae_vector vals;
ae_vector idx;
ae_vector ridx;
ae_vector didx;
ae_vector uidx;
ae_int_t matrixtype;
ae_int_t m;
ae_int_t n;
ae_int_t nfree;
ae_int_t ninitialized;
} sparsematrix;
typedef struct
{
ae_int_t n;
ae_int_t m;
ae_int_t nstart;
ae_int_t nits;
ae_int_t seedval;
ae_vector x0;
ae_vector x1;
ae_vector t;
ae_vector xbest;
hqrndstate r;
ae_vector x;
ae_vector mv;
ae_vector mtv;
ae_bool needmv;
ae_bool needmtv;
double repnorm;
rcommstate rstate;
} normestimatorstate;
} }
///////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////
// //
// THIS SECTION CONTAINS C++ INTERFACE // THIS SECTION CONTAINS C++ INTERFACE
// //
///////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////
namespace alglib namespace alglib
{ {
skipping to change at line 96 skipping to change at line 129
matinvreport(); matinvreport();
matinvreport(const matinvreport &rhs); matinvreport(const matinvreport &rhs);
matinvreport& operator=(const matinvreport &rhs); matinvreport& operator=(const matinvreport &rhs);
virtual ~matinvreport(); virtual ~matinvreport();
double &r1; double &r1;
double &rinf; double &rinf;
}; };
/************************************************************************* /*************************************************************************
Sparse matrix
You should use ALGLIB functions to work with sparse matrix.
Never try to access its fields directly!
*************************************************************************/
class _sparsematrix_owner
{
public:
_sparsematrix_owner();
_sparsematrix_owner(const _sparsematrix_owner &rhs);
_sparsematrix_owner& operator=(const _sparsematrix_owner &rhs);
virtual ~_sparsematrix_owner();
alglib_impl::sparsematrix* c_ptr();
alglib_impl::sparsematrix* c_ptr() const;
protected:
alglib_impl::sparsematrix *p_struct;
};
class sparsematrix : public _sparsematrix_owner
{
public:
sparsematrix();
sparsematrix(const sparsematrix &rhs);
sparsematrix& operator=(const sparsematrix &rhs);
virtual ~sparsematrix();
};
/*************************************************************************
This object stores state of the iterative norm estimation algorithm.
You should use ALGLIB functions to work with this object.
*************************************************************************/
class _normestimatorstate_owner
{
public:
_normestimatorstate_owner();
_normestimatorstate_owner(const _normestimatorstate_owner &rhs);
_normestimatorstate_owner& operator=(const _normestimatorstate_owner &r
hs);
virtual ~_normestimatorstate_owner();
alglib_impl::normestimatorstate* c_ptr();
alglib_impl::normestimatorstate* c_ptr() const;
protected:
alglib_impl::normestimatorstate *p_struct;
};
class normestimatorstate : public _normestimatorstate_owner
{
public:
normestimatorstate();
normestimatorstate(const normestimatorstate &rhs);
normestimatorstate& operator=(const normestimatorstate &rhs);
virtual ~normestimatorstate();
};
/*************************************************************************
Cache-oblivous complex "copy-and-transpose" Cache-oblivous complex "copy-and-transpose"
Input parameters: Input parameters:
M - number of rows M - number of rows
N - number of columns N - number of columns
A - source matrix, MxN submatrix is copied and transposed A - source matrix, MxN submatrix is copied and transposed
IA - submatrix offset (row index) IA - submatrix offset (row index)
JA - submatrix offset (column index) JA - submatrix offset (column index)
A - destination matrix B - destination matrix, must be large enough to store result
IB - submatrix offset (row index) IB - submatrix offset (row index)
JB - submatrix offset (column index) JB - submatrix offset (column index)
*************************************************************************/ *************************************************************************/
void cmatrixtranspose(const ae_int_t m, const ae_int_t n, const complex_2d_ array &a, const ae_int_t ia, const ae_int_t ja, complex_2d_array &b, const ae_int_t ib, const ae_int_t jb); void cmatrixtranspose(const ae_int_t m, const ae_int_t n, const complex_2d_ array &a, const ae_int_t ia, const ae_int_t ja, complex_2d_array &b, const ae_int_t ib, const ae_int_t jb);
/************************************************************************* /*************************************************************************
Cache-oblivous real "copy-and-transpose" Cache-oblivous real "copy-and-transpose"
Input parameters: Input parameters:
M - number of rows M - number of rows
N - number of columns N - number of columns
A - source matrix, MxN submatrix is copied and transposed A - source matrix, MxN submatrix is copied and transposed
IA - submatrix offset (row index) IA - submatrix offset (row index)
JA - submatrix offset (column index) JA - submatrix offset (column index)
A - destination matrix B - destination matrix, must be large enough to store result
IB - submatrix offset (row index) IB - submatrix offset (row index)
JB - submatrix offset (column index) JB - submatrix offset (column index)
*************************************************************************/ *************************************************************************/
void rmatrixtranspose(const ae_int_t m, const ae_int_t n, const real_2d_arr ay &a, const ae_int_t ia, const ae_int_t ja, real_2d_array &b, const ae_int _t ib, const ae_int_t jb); void rmatrixtranspose(const ae_int_t m, const ae_int_t n, const real_2d_arr ay &a, const ae_int_t ia, const ae_int_t ja, real_2d_array &b, const ae_int _t ib, const ae_int_t jb);
/************************************************************************* /*************************************************************************
Copy Copy
Input parameters: Input parameters:
M - number of rows M - number of rows
N - number of columns N - number of columns
A - source matrix, MxN submatrix is copied and transposed A - source matrix, MxN submatrix is copied and transposed
IA - submatrix offset (row index) IA - submatrix offset (row index)
JA - submatrix offset (column index) JA - submatrix offset (column index)
B - destination matrix B - destination matrix, must be large enough to store result
IB - submatrix offset (row index) IB - submatrix offset (row index)
JB - submatrix offset (column index) JB - submatrix offset (column index)
*************************************************************************/ *************************************************************************/
void cmatrixcopy(const ae_int_t m, const ae_int_t n, const complex_2d_array &a, const ae_int_t ia, const ae_int_t ja, complex_2d_array &b, const ae_in t_t ib, const ae_int_t jb); void cmatrixcopy(const ae_int_t m, const ae_int_t n, const complex_2d_array &a, const ae_int_t ia, const ae_int_t ja, complex_2d_array &b, const ae_in t_t ib, const ae_int_t jb);
/************************************************************************* /*************************************************************************
Copy Copy
Input parameters: Input parameters:
M - number of rows M - number of rows
N - number of columns N - number of columns
A - source matrix, MxN submatrix is copied and transposed A - source matrix, MxN submatrix is copied and transposed
IA - submatrix offset (row index) IA - submatrix offset (row index)
JA - submatrix offset (column index) JA - submatrix offset (column index)
B - destination matrix B - destination matrix, must be large enough to store result
IB - submatrix offset (row index) IB - submatrix offset (row index)
JB - submatrix offset (column index) JB - submatrix offset (column index)
*************************************************************************/ *************************************************************************/
void rmatrixcopy(const ae_int_t m, const ae_int_t n, const real_2d_array &a , const ae_int_t ia, const ae_int_t ja, real_2d_array &b, const ae_int_t ib , const ae_int_t jb); void rmatrixcopy(const ae_int_t m, const ae_int_t n, const real_2d_array &a , const ae_int_t ia, const ae_int_t ja, real_2d_array &b, const ae_int_t ib , const ae_int_t jb);
/************************************************************************* /*************************************************************************
Rank-1 correction: A := A + u*v' Rank-1 correction: A := A + u*v'
INPUT PARAMETERS: INPUT PARAMETERS:
M - number of rows M - number of rows
skipping to change at line 205 skipping to change at line 293
A - target matrix A - target matrix
IA - submatrix offset (row index) IA - submatrix offset (row index)
JA - submatrix offset (column index) JA - submatrix offset (column index)
OpA - operation type: OpA - operation type:
* OpA=0 => op(A) = A * OpA=0 => op(A) = A
* OpA=1 => op(A) = A^T * OpA=1 => op(A) = A^T
* OpA=2 => op(A) = A^H * OpA=2 => op(A) = A^H
X - input vector X - input vector
IX - subvector offset IX - subvector offset
IY - subvector offset IY - subvector offset
Y - preallocated matrix, must be large enough to store result
OUTPUT PARAMETERS: OUTPUT PARAMETERS:
Y - vector which stores result Y - vector which stores result
if M=0, then subroutine does nothing. if M=0, then subroutine does nothing.
if N=0, Y is filled by zeros. if N=0, Y is filled by zeros.
-- ALGLIB routine -- -- ALGLIB routine --
28.01.2010 28.01.2010
skipping to change at line 234 skipping to change at line 323
N - number of columns of op(A) N - number of columns of op(A)
A - target matrix A - target matrix
IA - submatrix offset (row index) IA - submatrix offset (row index)
JA - submatrix offset (column index) JA - submatrix offset (column index)
OpA - operation type: OpA - operation type:
* OpA=0 => op(A) = A * OpA=0 => op(A) = A
* OpA=1 => op(A) = A^T * OpA=1 => op(A) = A^T
X - input vector X - input vector
IX - subvector offset IX - subvector offset
IY - subvector offset IY - subvector offset
Y - preallocated matrix, must be large enough to store result
OUTPUT PARAMETERS: OUTPUT PARAMETERS:
Y - vector which stores result Y - vector which stores result
if M=0, then subroutine does nothing. if M=0, then subroutine does nothing.
if N=0, Y is filled by zeros. if N=0, Y is filled by zeros.
-- ALGLIB routine -- -- ALGLIB routine --
28.01.2010 28.01.2010
skipping to change at line 269 skipping to change at line 359
M - matrix size, N>=0 M - matrix size, N>=0
A - matrix, actial matrix is stored in A[I1:I1+N-1,J1:J1+N-1] A - matrix, actial matrix is stored in A[I1:I1+N-1,J1:J1+N-1]
I1 - submatrix offset I1 - submatrix offset
J1 - submatrix offset J1 - submatrix offset
IsUpper - whether matrix is upper triangular IsUpper - whether matrix is upper triangular
IsUnit - whether matrix is unitriangular IsUnit - whether matrix is unitriangular
OpType - transformation type: OpType - transformation type:
* 0 - no transformation * 0 - no transformation
* 1 - transposition * 1 - transposition
* 2 - conjugate transposition * 2 - conjugate transposition
C - matrix, actial matrix is stored in C[I2:I2+M-1,J2:J2+N-1] X - matrix, actial matrix is stored in X[I2:I2+M-1,J2:J2+N-1]
I2 - submatrix offset I2 - submatrix offset
J2 - submatrix offset J2 - submatrix offset
-- ALGLIB routine -- -- ALGLIB routine --
15.12.2009 15.12.2009
Bochkanov Sergey Bochkanov Sergey
*************************************************************************/ *************************************************************************/
void cmatrixrighttrsm(const ae_int_t m, const ae_int_t n, const complex_2d_ array &a, const ae_int_t i1, const ae_int_t j1, const bool isupper, const b ool isunit, const ae_int_t optype, complex_2d_array &x, const ae_int_t i2, const ae_int_t j2); void cmatrixrighttrsm(const ae_int_t m, const ae_int_t n, const complex_2d_ array &a, const ae_int_t i1, const ae_int_t j1, const bool isupper, const b ool isunit, const ae_int_t optype, complex_2d_array &x, const ae_int_t i2, const ae_int_t j2);
/************************************************************************* /*************************************************************************
skipping to change at line 300 skipping to change at line 390
M - matrix size, N>=0 M - matrix size, N>=0
A - matrix, actial matrix is stored in A[I1:I1+M-1,J1:J1+M-1] A - matrix, actial matrix is stored in A[I1:I1+M-1,J1:J1+M-1]
I1 - submatrix offset I1 - submatrix offset
J1 - submatrix offset J1 - submatrix offset
IsUpper - whether matrix is upper triangular IsUpper - whether matrix is upper triangular
IsUnit - whether matrix is unitriangular IsUnit - whether matrix is unitriangular
OpType - transformation type: OpType - transformation type:
* 0 - no transformation * 0 - no transformation
* 1 - transposition * 1 - transposition
* 2 - conjugate transposition * 2 - conjugate transposition
C - matrix, actial matrix is stored in C[I2:I2+M-1,J2:J2+N-1] X - matrix, actial matrix is stored in X[I2:I2+M-1,J2:J2+N-1]
I2 - submatrix offset I2 - submatrix offset
J2 - submatrix offset J2 - submatrix offset
-- ALGLIB routine -- -- ALGLIB routine --
15.12.2009 15.12.2009
Bochkanov Sergey Bochkanov Sergey
*************************************************************************/ *************************************************************************/
void cmatrixlefttrsm(const ae_int_t m, const ae_int_t n, const complex_2d_a rray &a, const ae_int_t i1, const ae_int_t j1, const bool isupper, const bo ol isunit, const ae_int_t optype, complex_2d_array &x, const ae_int_t i2, c onst ae_int_t j2); void cmatrixlefttrsm(const ae_int_t m, const ae_int_t n, const complex_2d_a rray &a, const ae_int_t i1, const ae_int_t j1, const bool isupper, const bo ol isunit, const ae_int_t optype, complex_2d_array &x, const ae_int_t i2, c onst ae_int_t j2);
/************************************************************************* /*************************************************************************
Same as CMatrixRightTRSM, but for real matrices This subroutine calculates X*op(A^-1) where:
* X is MxN general matrix
* A is NxN upper/lower triangular/unitriangular matrix
* "op" may be identity transformation, transposition
Multiplication result replaces X.
Cache-oblivious algorithm is used.
OpType may be only 0 or 1. INPUT PARAMETERS
N - matrix size, N>=0
M - matrix size, N>=0
A - matrix, actial matrix is stored in A[I1:I1+N-1,J1:J1+N-1]
I1 - submatrix offset
J1 - submatrix offset
IsUpper - whether matrix is upper triangular
IsUnit - whether matrix is unitriangular
OpType - transformation type:
* 0 - no transformation
* 1 - transposition
X - matrix, actial matrix is stored in X[I2:I2+M-1,J2:J2+N-1]
I2 - submatrix offset
J2 - submatrix offset
-- ALGLIB routine -- -- ALGLIB routine --
15.12.2009 15.12.2009
Bochkanov Sergey Bochkanov Sergey
*************************************************************************/ *************************************************************************/
void rmatrixrighttrsm(const ae_int_t m, const ae_int_t n, const real_2d_arr ay &a, const ae_int_t i1, const ae_int_t j1, const bool isupper, const bool isunit, const ae_int_t optype, real_2d_array &x, const ae_int_t i2, const ae_int_t j2); void rmatrixrighttrsm(const ae_int_t m, const ae_int_t n, const real_2d_arr ay &a, const ae_int_t i1, const ae_int_t j1, const bool isupper, const bool isunit, const ae_int_t optype, real_2d_array &x, const ae_int_t i2, const ae_int_t j2);
/************************************************************************* /*************************************************************************
Same as CMatrixLeftTRSM, but for real matrices This subroutine calculates op(A^-1)*X where:
* X is MxN general matrix
* A is MxM upper/lower triangular/unitriangular matrix
* "op" may be identity transformation, transposition
Multiplication result replaces X.
Cache-oblivious algorithm is used.
OpType may be only 0 or 1. INPUT PARAMETERS
N - matrix size, N>=0
M - matrix size, N>=0
A - matrix, actial matrix is stored in A[I1:I1+M-1,J1:J1+M-1]
I1 - submatrix offset
J1 - submatrix offset
IsUpper - whether matrix is upper triangular
IsUnit - whether matrix is unitriangular
OpType - transformation type:
* 0 - no transformation
* 1 - transposition
X - matrix, actial matrix is stored in X[I2:I2+M-1,J2:J2+N-1]
I2 - submatrix offset
J2 - submatrix offset
-- ALGLIB routine -- -- ALGLIB routine --
15.12.2009 15.12.2009
Bochkanov Sergey Bochkanov Sergey
*************************************************************************/ *************************************************************************/
void rmatrixlefttrsm(const ae_int_t m, const ae_int_t n, const real_2d_arra y &a, const ae_int_t i1, const ae_int_t j1, const bool isupper, const bool isunit, const ae_int_t optype, real_2d_array &x, const ae_int_t i2, const a e_int_t j2); void rmatrixlefttrsm(const ae_int_t m, const ae_int_t n, const real_2d_arra y &a, const ae_int_t i1, const ae_int_t j1, const bool isupper, const bool isunit, const ae_int_t optype, real_2d_array &x, const ae_int_t i2, const a e_int_t j2);
/************************************************************************* /*************************************************************************
This subroutine calculates C=alpha*A*A^H+beta*C or C=alpha*A^H*A+beta*C This subroutine calculates C=alpha*A*A^H+beta*C or C=alpha*A^H*A+beta*C
where: where:
skipping to change at line 368 skipping to change at line 496
JC - submatrix offset JC - submatrix offset
IsUpper - whether C is upper triangular or lower triangular IsUpper - whether C is upper triangular or lower triangular
-- ALGLIB routine -- -- ALGLIB routine --
16.12.2009 16.12.2009
Bochkanov Sergey Bochkanov Sergey
*************************************************************************/ *************************************************************************/
void cmatrixsyrk(const ae_int_t n, const ae_int_t k, const double alpha, co nst complex_2d_array &a, const ae_int_t ia, const ae_int_t ja, const ae_int _t optypea, const double beta, complex_2d_array &c, const ae_int_t ic, cons t ae_int_t jc, const bool isupper); void cmatrixsyrk(const ae_int_t n, const ae_int_t k, const double alpha, co nst complex_2d_array &a, const ae_int_t ia, const ae_int_t ja, const ae_int _t optypea, const double beta, complex_2d_array &c, const ae_int_t ic, cons t ae_int_t jc, const bool isupper);
/************************************************************************* /*************************************************************************
Same as CMatrixSYRK, but for real matrices This subroutine calculates C=alpha*A*A^T+beta*C or C=alpha*A^T*A+beta*C
where:
* C is NxN symmetric matrix given by its upper/lower triangle
* A is NxK matrix when A*A^T is calculated, KxN matrix otherwise
Additional info:
* cache-oblivious algorithm is used.
* multiplication result replaces C. If Beta=0, C elements are not used in
calculations (not multiplied by zero - just not referenced)
* if Alpha=0, A is not used (not multiplied by zero - just not referenced)
* if both Beta and Alpha are zero, C is filled by zeros.
OpType may be only 0 or 1. INPUT PARAMETERS
N - matrix size, N>=0
K - matrix size, K>=0
Alpha - coefficient
A - matrix
IA - submatrix offset
JA - submatrix offset
OpTypeA - multiplication type:
* 0 - A*A^T is calculated
* 2 - A^T*A is calculated
Beta - coefficient
C - matrix
IC - submatrix offset
JC - submatrix offset
IsUpper - whether C is upper triangular or lower triangular
-- ALGLIB routine -- -- ALGLIB routine --
16.12.2009 16.12.2009
Bochkanov Sergey Bochkanov Sergey
*************************************************************************/ *************************************************************************/
void rmatrixsyrk(const ae_int_t n, const ae_int_t k, const double alpha, co nst real_2d_array &a, const ae_int_t ia, const ae_int_t ja, const ae_int_t optypea, const double beta, real_2d_array &c, const ae_int_t ic, const ae_i nt_t jc, const bool isupper); void rmatrixsyrk(const ae_int_t n, const ae_int_t k, const double alpha, co nst real_2d_array &a, const ae_int_t ia, const ae_int_t ja, const ae_int_t optypea, const double beta, real_2d_array &c, const ae_int_t ic, const ae_i nt_t jc, const bool isupper);
/************************************************************************* /*************************************************************************
This subroutine calculates C = alpha*op1(A)*op2(B) +beta*C where: This subroutine calculates C = alpha*op1(A)*op2(B) +beta*C where:
* C is MxN general matrix * C is MxN general matrix
skipping to change at line 393 skipping to change at line 545
* "op" may be identity transformation, transposition, conjugate transpositi on * "op" may be identity transformation, transposition, conjugate transpositi on
Additional info: Additional info:
* cache-oblivious algorithm is used. * cache-oblivious algorithm is used.
* multiplication result replaces C. If Beta=0, C elements are not used in * multiplication result replaces C. If Beta=0, C elements are not used in
calculations (not multiplied by zero - just not referenced) calculations (not multiplied by zero - just not referenced)
* if Alpha=0, A is not used (not multiplied by zero - just not referenced) * if Alpha=0, A is not used (not multiplied by zero - just not referenced)
* if both Beta and Alpha are zero, C is filled by zeros. * if both Beta and Alpha are zero, C is filled by zeros.
INPUT PARAMETERS INPUT PARAMETERS
M - matrix size, M>0
N - matrix size, N>0 N - matrix size, N>0
M - matrix size, N>0
K - matrix size, K>0 K - matrix size, K>0
Alpha - coefficient Alpha - coefficient
A - matrix A - matrix
IA - submatrix offset IA - submatrix offset
JA - submatrix offset JA - submatrix offset
OpTypeA - transformation type: OpTypeA - transformation type:
* 0 - no transformation * 0 - no transformation
* 1 - transposition * 1 - transposition
* 2 - conjugate transposition * 2 - conjugate transposition
B - matrix B - matrix
skipping to change at line 423 skipping to change at line 575
IC - submatrix offset IC - submatrix offset
JC - submatrix offset JC - submatrix offset
-- ALGLIB routine -- -- ALGLIB routine --
16.12.2009 16.12.2009
Bochkanov Sergey Bochkanov Sergey
*************************************************************************/ *************************************************************************/
void cmatrixgemm(const ae_int_t m, const ae_int_t n, const ae_int_t k, cons t alglib::complex alpha, const complex_2d_array &a, const ae_int_t ia, cons t ae_int_t ja, const ae_int_t optypea, const complex_2d_array &b, const ae_ int_t ib, const ae_int_t jb, const ae_int_t optypeb, const alglib::complex beta, complex_2d_array &c, const ae_int_t ic, const ae_int_t jc); void cmatrixgemm(const ae_int_t m, const ae_int_t n, const ae_int_t k, cons t alglib::complex alpha, const complex_2d_array &a, const ae_int_t ia, cons t ae_int_t ja, const ae_int_t optypea, const complex_2d_array &b, const ae_ int_t ib, const ae_int_t jb, const ae_int_t optypeb, const alglib::complex beta, complex_2d_array &c, const ae_int_t ic, const ae_int_t jc);
/************************************************************************* /*************************************************************************
Same as CMatrixGEMM, but for real numbers. This subroutine calculates C = alpha*op1(A)*op2(B) +beta*C where:
OpType may be only 0 or 1. * C is MxN general matrix
* op1(A) is MxK matrix
* op2(B) is KxN matrix
* "op" may be identity transformation, transposition
Additional info:
* cache-oblivious algorithm is used.
* multiplication result replaces C. If Beta=0, C elements are not used in
calculations (not multiplied by zero - just not referenced)
* if Alpha=0, A is not used (not multiplied by zero - just not referenced)
* if both Beta and Alpha are zero, C is filled by zeros.
INPUT PARAMETERS
M - matrix size, M>0
N - matrix size, N>0
K - matrix size, K>0
Alpha - coefficient
A - matrix
IA - submatrix offset
JA - submatrix offset
OpTypeA - transformation type:
* 0 - no transformation
* 1 - transposition
B - matrix
IB - submatrix offset
JB - submatrix offset
OpTypeB - transformation type:
* 0 - no transformation
* 1 - transposition
Beta - coefficient
C - matrix
IC - submatrix offset
JC - submatrix offset
-- ALGLIB routine -- -- ALGLIB routine --
16.12.2009 16.12.2009
Bochkanov Sergey Bochkanov Sergey
*************************************************************************/ *************************************************************************/
void rmatrixgemm(const ae_int_t m, const ae_int_t n, const ae_int_t k, cons t double alpha, const real_2d_array &a, const ae_int_t ia, const ae_int_t j a, const ae_int_t optypea, const real_2d_array &b, const ae_int_t ib, const ae_int_t jb, const ae_int_t optypeb, const double beta, real_2d_array &c, const ae_int_t ic, const ae_int_t jc); void rmatrixgemm(const ae_int_t m, const ae_int_t n, const ae_int_t k, cons t double alpha, const real_2d_array &a, const ae_int_t ia, const ae_int_t j a, const ae_int_t optypea, const real_2d_array &b, const ae_int_t ib, const ae_int_t jb, const ae_int_t optypeb, const double beta, real_2d_array &c, const ae_int_t ic, const ae_int_t jc);
/************************************************************************* /*************************************************************************
QR decomposition of a rectangular matrix of size MxN QR decomposition of a rectangular matrix of size MxN
skipping to change at line 1162 skipping to change at line 1346
Output parameters: Output parameters:
Q - transformation matrix. Q - transformation matrix.
array with elements [0..N-1, 0..N-1]. array with elements [0..N-1, 0..N-1].
-- ALGLIB -- -- ALGLIB --
Copyright 2005-2010 by Bochkanov Sergey Copyright 2005-2010 by Bochkanov Sergey
*************************************************************************/ *************************************************************************/
void hmatrixtdunpackq(const complex_2d_array &a, const ae_int_t n, const bo ol isupper, const complex_1d_array &tau, complex_2d_array &q); void hmatrixtdunpackq(const complex_2d_array &a, const ae_int_t n, const bo ol isupper, const complex_1d_array &tau, complex_2d_array &q);
/************************************************************************* /*************************************************************************
Singular value decomposition of a bidiagonal matrix (extended algorithm)
The algorithm performs the singular value decomposition of a bidiagonal
matrix B (upper or lower) representing it as B = Q*S*P^T, where Q and P -
orthogonal matrices, S - diagonal matrix with non-negative elements on the
main diagonal, in descending order.
The algorithm finds singular values. In addition, the algorithm can
calculate matrices Q and P (more precisely, not the matrices, but their
product with given matrices U and VT - U*Q and (P^T)*VT)). Of course,
matrices U and VT can be of any type, including identity. Furthermore, the
algorithm can calculate Q'*C (this product is calculated more effectively
than U*Q, because this calculation operates with rows instead of matrix
columns).
The feature of the algorithm is its ability to find all singular values
including those which are arbitrarily close to 0 with relative accuracy
close to machine precision. If the parameter IsFractionalAccuracyRequired
is set to True, all singular values will have high relative accuracy close
to machine precision. If the parameter is set to False, only the biggest
singular value will have relative accuracy close to machine precision.
The absolute error of other singular values is equal to the absolute error
of the biggest singular value.
Input parameters:
D - main diagonal of matrix B.
Array whose index ranges within [0..N-1].
E - superdiagonal (or subdiagonal) of matrix B.
Array whose index ranges within [0..N-2].
N - size of matrix B.
IsUpper - True, if the matrix is upper bidiagonal.
IsFractionalAccuracyRequired -
THIS PARAMETER IS IGNORED SINCE ALGLIB 3.5.0
SINGULAR VALUES ARE ALWAYS SEARCHED WITH HIGH ACCURACY.
U - matrix to be multiplied by Q.
Array whose indexes range within [0..NRU-1, 0..N-1].
The matrix can be bigger, in that case only the submatrix
[0..NRU-1, 0..N-1] will be multiplied by Q.
NRU - number of rows in matrix U.
C - matrix to be multiplied by Q'.
Array whose indexes range within [0..N-1, 0..NCC-1].
The matrix can be bigger, in that case only the submatrix
[0..N-1, 0..NCC-1] will be multiplied by Q'.
NCC - number of columns in matrix C.
VT - matrix to be multiplied by P^T.
Array whose indexes range within [0..N-1, 0..NCVT-1].
The matrix can be bigger, in that case only the submatrix
[0..N-1, 0..NCVT-1] will be multiplied by P^T.
NCVT - number of columns in matrix VT.
Output parameters:
D - singular values of matrix B in descending order.
U - if NRU>0, contains matrix U*Q.
VT - if NCVT>0, contains matrix (P^T)*VT.
C - if NCC>0, contains matrix Q'*C.
Result:
True, if the algorithm has converged.
False, if the algorithm hasn't converged (rare case).
Additional information:
The type of convergence is controlled by the internal parameter TOL.
If the parameter is greater than 0, the singular values will have
relative accuracy TOL. If TOL<0, the singular values will have
absolute accuracy ABS(TOL)*norm(B).
By default, |TOL| falls within the range of 10*Epsilon and 100*Epsilon,
where Epsilon is the machine precision. It is not recommended to use
TOL less than 10*Epsilon since this will considerably slow down the
algorithm and may not lead to error decreasing.
History:
* 31 March, 2007.
changed MAXITR from 6 to 12.
-- LAPACK routine (version 3.0) --
Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
Courant Institute, Argonne National Lab, and Rice University
October 31, 1999.
*************************************************************************/
bool rmatrixbdsvd(real_1d_array &d, const real_1d_array &e, const ae_int_t
n, const bool isupper, const bool isfractionalaccuracyrequired, real_2d_arr
ay &u, const ae_int_t nru, real_2d_array &c, const ae_int_t ncc, real_2d_ar
ray &vt, const ae_int_t ncvt);
/*************************************************************************
Singular value decomposition of a rectangular matrix.
The algorithm calculates the singular value decomposition of a matrix of
size MxN: A = U * S * V^T
The algorithm finds the singular values and, optionally, matrices U and V^T
.
The algorithm can find both first min(M,N) columns of matrix U and rows of
matrix V^T (singular vectors), and matrices U and V^T wholly (of sizes MxM
and NxN respectively).
Take into account that the subroutine does not return matrix V but V^T.
Input parameters:
A - matrix to be decomposed.
Array whose indexes range within [0..M-1, 0..N-1].
M - number of rows in matrix A.
N - number of columns in matrix A.
UNeeded - 0, 1 or 2. See the description of the parameter U.
VTNeeded - 0, 1 or 2. See the description of the parameter VT.
AdditionalMemory -
If the parameter:
* equals 0, the algorithm doesn
memory (lower requirements, lower performance).
* equals 1, the algorithm uses additional
memory of size min(M,N)*min(M,N) of real numbers.
It often speeds up the algorithm.
* equals 2, the algorithm uses additional
memory of size M*min(M,N) of real numbers.
It allows to get a maximum performance.
The recommended value of the parameter is 2.
Output parameters:
W - contains singular values in descending order.
U - if UNeeded=0, U isn't changed, the left singular vector
s
are not calculated.
if Uneeded=1, U contains left singular vectors (first
min(M,N) columns of matrix U). Array whose indexes rang
e
within [0..M-1, 0..Min(M,N)-1].
if UNeeded=2, U contains matrix U wholly. Array whose
indexes range within [0..M-1, 0..M-1].
VT - if VTNeeded=0, VT isn
are not calculated.
if VTNeeded=1, VT contains right singular vectors (firs
t
min(M,N) rows of matrix V^T). Array whose indexes range
within [0..min(M,N)-1, 0..N-1].
if VTNeeded=2, VT contains matrix V^T wholly. Array who
se
indexes range within [0..N-1, 0..N-1].
-- ALGLIB --
Copyright 2005 by Bochkanov Sergey
*************************************************************************/
bool rmatrixsvd(const real_2d_array &a, const ae_int_t m, const ae_int_t n,
const ae_int_t uneeded, const ae_int_t vtneeded, const ae_int_t additional
memory, real_1d_array &w, real_2d_array &u, real_2d_array &vt);
/*************************************************************************
Finding the eigenvalues and eigenvectors of a symmetric matrix Finding the eigenvalues and eigenvectors of a symmetric matrix
The algorithm finds eigen pairs of a symmetric matrix by reducing it to The algorithm finds eigen pairs of a symmetric matrix by reducing it to
tridiagonal form and using the QL/QR algorithm. tridiagonal form and using the QL/QR algorithm.
Input parameters: Input parameters:
A - symmetric matrix which is given by its upper or lower A - symmetric matrix which is given by its upper or lower
triangular part. triangular part.
Array whose indexes range within [0..N-1, 0..N-1]. Array whose indexes range within [0..N-1, 0..N-1].
N - size of matrix A. N - size of matrix A.
skipping to change at line 2679 skipping to change at line 2998
Rep - same as for RMatrixLUInverse Rep - same as for RMatrixLUInverse
A - same as for RMatrixLUInverse. A - same as for RMatrixLUInverse.
-- ALGLIB -- -- ALGLIB --
Copyright 05.02.2010 by Bochkanov Sergey Copyright 05.02.2010 by Bochkanov Sergey
*************************************************************************/ *************************************************************************/
void cmatrixtrinverse(complex_2d_array &a, const ae_int_t n, const bool isu pper, const bool isunit, ae_int_t &info, matinvreport &rep); void cmatrixtrinverse(complex_2d_array &a, const ae_int_t n, const bool isu pper, const bool isunit, ae_int_t &info, matinvreport &rep);
void cmatrixtrinverse(complex_2d_array &a, const bool isupper, ae_int_t &in fo, matinvreport &rep); void cmatrixtrinverse(complex_2d_array &a, const bool isupper, ae_int_t &in fo, matinvreport &rep);
/************************************************************************* /*************************************************************************
Singular value decomposition of a bidiagonal matrix (extended algorithm) This function creates sparse matrix in a Hash-Table format.
The algorithm performs the singular value decomposition of a bidiagonal This function creates Hast-Table matrix, which can be converted to CRS
matrix B (upper or lower) representing it as B = Q*S*P^T, where Q and P - format after its initialization is over. Typical usage scenario for a
orthogonal matrices, S - diagonal matrix with non-negative elements on the sparse matrix is:
main diagonal, in descending order. 1. creation in a Hash-Table format
2. insertion of the matrix elements
3. conversion to the CRS representation
4. matrix is passed to some linear algebra algorithm
The algorithm finds singular values. In addition, the algorithm can Some information about different matrix formats can be found below, in
calculate matrices Q and P (more precisely, not the matrices, but their the "NOTES" section.
product with given matrices U and VT - U*Q and (P^T)*VT)). Of course,
matrices U and VT can be of any type, including identity. Furthermore, the
algorithm can calculate Q'*C (this product is calculated more effectively
than U*Q, because this calculation operates with rows instead of matrix
columns).
The feature of the algorithm is its ability to find all singular values INPUT PARAMETERS
including those which are arbitrarily close to 0 with relative accuracy M - number of rows in a matrix, M>=1
close to machine precision. If the parameter IsFractionalAccuracyRequired N - number of columns in a matrix, N>=1
is set to True, all singular values will have high relative accuracy close K - K>=0, expected number of non-zero elements in a matrix.
to machine precision. If the parameter is set to False, only the biggest K can be inexact approximation, can be less than actual
singular value will have relative accuracy close to machine precision. number of elements (table will grow when needed) or
The absolute error of other singular values is equal to the absolute error even zero).
of the biggest singular value. It is important to understand that although hash-table
may grow automatically, it is better to provide good
estimate of data size.
OUTPUT PARAMETERS
S - sparse M*N matrix in Hash-Table representation.
All elements of the matrix are zero.
NOTE 1.
Sparse matrices can be stored using either Hash-Table representation or
Compressed Row Storage representation. Hast-table is better suited for
querying and dynamic operations (thus, it is used for matrix
initialization), but it is inefficient when you want to make some linear
algebra operations.
From the other side, CRS is better suited for linear algebra operations,
but initialization is less convenient - you have to tell row sizes at the
initialization, and you can fill matrix only row by row, from left to
right. CRS is also very inefficient when you want to find matrix element
by its index.
Thus, Hash-Table representation does not support linear algebra
operations, while CRS format does not support modification of the table.
Tables below outline information about these two formats:
OPERATIONS WITH MATRIX HASH CRS
create + +
read element + +
modify element +
add value to element +
A*x (dense vector) +
A'*x (dense vector) +
A*X (dense matrix) +
A'*X (dense matrix) +
NOTE 2.
Hash-tables use memory inefficiently, and they have to keep some amount
of the "spare memory" in order to have good performance. Hash table for
matrix with K non-zero elements will need C*K*(8+2*sizeof(int)) bytes,
where C is a small constant, about 1.5-2 in magnitude.
CRS storage, from the other side, is more memory-efficient, and needs
just K*(8+sizeof(int))+M*sizeof(int) bytes, where M is a number of rows
in a matrix.
When you convert from the Hash-Table to CRS representation, all unneeded
memory will be freed.
-- ALGLIB PROJECT --
Copyright 14.10.2011 by Bochkanov Sergey
*************************************************************************/
void sparsecreate(const ae_int_t m, const ae_int_t n, const ae_int_t k, spa
rsematrix &s);
void sparsecreate(const ae_int_t m, const ae_int_t n, sparsematrix &s);
/*************************************************************************
This function creates sparse matrix in a CRS format (expert function for
situations when you are running out of memory).
This function creates CRS matrix. Typical usage scenario for a CRS matrix
is:
1. creation (you have to tell number of non-zero elements at each row at
this moment)
2. insertion of the matrix elements (row by row, from left to right)
3. matrix is passed to some linear algebra algorithm
This function is a memory-efficient alternative to SparseCreate(), but it
is more complex because it requires you to know in advance how large your
matrix is. Some information about different matrix formats can be found
below, in the "NOTES" section.
Input parameters: INPUT PARAMETERS
D - main diagonal of matrix B. M - number of rows in a matrix, M>=1
Array whose index ranges within [0..N-1]. N - number of columns in a matrix, N>=1
E - superdiagonal (or subdiagonal) of matrix B. NER - number of elements at each row, array[M], NER[i]>=0
Array whose index ranges within [0..N-2].
N - size of matrix B.
IsUpper - True, if the matrix is upper bidiagonal.
IsFractionalAccuracyRequired -
accuracy to search singular values with.
U - matrix to be multiplied by Q.
Array whose indexes range within [0..NRU-1, 0..N-1].
The matrix can be bigger, in that case only the submatrix
[0..NRU-1, 0..N-1] will be multiplied by Q.
NRU - number of rows in matrix U.
C - matrix to be multiplied by Q'.
Array whose indexes range within [0..N-1, 0..NCC-1].
The matrix can be bigger, in that case only the submatrix
[0..N-1, 0..NCC-1] will be multiplied by Q'.
NCC - number of columns in matrix C.
VT - matrix to be multiplied by P^T.
Array whose indexes range within [0..N-1, 0..NCVT-1].
The matrix can be bigger, in that case only the submatrix
[0..N-1, 0..NCVT-1] will be multiplied by P^T.
NCVT - number of columns in matrix VT.
Output parameters: OUTPUT PARAMETERS
D - singular values of matrix B in descending order. S - sparse M*N matrix in CRS representation.
U - if NRU>0, contains matrix U*Q. You have to fill ALL non-zero elements by calling
VT - if NCVT>0, contains matrix (P^T)*VT. SparseSet() BEFORE you try to use this matrix.
C - if NCC>0, contains matrix Q'*C.
Result: NOTE 1.
True, if the algorithm has converged.
False, if the algorithm hasn't converged (rare case).
Additional information: Sparse matrices can be stored using either Hash-Table representation or
The type of convergence is controlled by the internal parameter TOL. Compressed Row Storage representation. Hast-table is better suited for
If the parameter is greater than 0, the singular values will have querying and dynamic operations (thus, it is used for matrix
relative accuracy TOL. If TOL<0, the singular values will have initialization), but it is inefficient when you want to make some linear
absolute accuracy ABS(TOL)*norm(B). algebra operations.
By default, |TOL| falls within the range of 10*Epsilon and 100*Epsilon,
where Epsilon is the machine precision. It is not recommended to use
TOL less than 10*Epsilon since this will considerably slow down the
algorithm and may not lead to error decreasing.
History:
* 31 March, 2007.
changed MAXITR from 6 to 12.
-- LAPACK routine (version 3.0) -- From the other side, CRS is better suited for linear algebra operations,
Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., but initialization is less convenient - you have to tell row sizes at the
Courant Institute, Argonne National Lab, and Rice University initialization, and you can fill matrix only row by row, from left to
October 31, 1999. right. CRS is also very inefficient when you want to find matrix element
by its index.
Thus, Hash-Table representation does not support linear algebra
operations, while CRS format does not support modification of the table.
Tables below outline information about these two formats:
OPERATIONS WITH MATRIX HASH CRS
create + +
read element + +
modify element +
add value to element +
A*x (dense vector) +
A'*x (dense vector) +
A*X (dense matrix) +
A'*X (dense matrix) +
NOTE 2.
Hash-tables use memory inefficiently, and they have to keep some amount
of the "spare memory" in order to have good performance. Hash table for
matrix with K non-zero elements will need C*K*(8+2*sizeof(int)) bytes,
where C is a small constant, about 1.5-2 in magnitude.
CRS storage, from the other side, is more memory-efficient, and needs
just K*(8+sizeof(int))+M*sizeof(int) bytes, where M is a number of rows
in a matrix.
When you convert from the Hash-Table to CRS representation, all unneeded
memory will be freed.
-- ALGLIB PROJECT --
Copyright 14.10.2011 by Bochkanov Sergey
*************************************************************************/ *************************************************************************/
bool rmatrixbdsvd(real_1d_array &d, const real_1d_array &e, const ae_int_t n, const bool isupper, const bool isfractionalaccuracyrequired, real_2d_arr ay &u, const ae_int_t nru, real_2d_array &c, const ae_int_t ncc, real_2d_ar ray &vt, const ae_int_t ncvt); void sparsecreatecrs(const ae_int_t m, const ae_int_t n, const integer_1d_a rray &ner, sparsematrix &s);
/************************************************************************* /*************************************************************************
Singular value decomposition of a rectangular matrix. This function copies S0 to S1.
The algorithm calculates the singular value decomposition of a matrix of NOTE: this function does not verify its arguments, it just copies all
size MxN: A = U * S * V^T fields of the structure.
The algorithm finds the singular values and, optionally, matrices U and V^T -- ALGLIB PROJECT --
. Copyright 14.10.2011 by Bochkanov Sergey
The algorithm can find both first min(M,N) columns of matrix U and rows of *************************************************************************/
matrix V^T (singular vectors), and matrices U and V^T wholly (of sizes MxM void sparsecopy(const sparsematrix &s0, sparsematrix &s1);
and NxN respectively).
Take into account that the subroutine does not return matrix V but V^T. /*************************************************************************
This function adds value to S[i,j] - element of the sparse matrix. Matrix
must be in a Hash-Table mode.
Input parameters: In case S[i,j] already exists in the table, V i added to its value. In
A - matrix to be decomposed. case S[i,j] is non-existent, it is inserted in the table. Table
Array whose indexes range within [0..M-1, 0..N-1]. automatically grows when necessary.
M - number of rows in matrix A.
N - number of columns in matrix A.
UNeeded - 0, 1 or 2. See the description of the parameter U.
VTNeeded - 0, 1 or 2. See the description of the parameter VT.
AdditionalMemory -
If the parameter:
* equals 0, the algorithm doesn
memory (lower requirements, lower performance).
* equals 1, the algorithm uses additional
memory of size min(M,N)*min(M,N) of real numbers.
It often speeds up the algorithm.
* equals 2, the algorithm uses additional
memory of size M*min(M,N) of real numbers.
It allows to get a maximum performance.
The recommended value of the parameter is 2.
Output parameters: INPUT PARAMETERS
W - contains singular values in descending order. S - sparse M*N matrix in Hash-Table representation.
U - if UNeeded=0, U isn't changed, the left singular vector Exception will be thrown for CRS matrix.
s I - row index of the element to modify, 0<=I<M
are not calculated. J - column index of the element to modify, 0<=J<N
if Uneeded=1, U contains left singular vectors (first V - value to add, must be finite number
min(M,N) columns of matrix U). Array whose indexes rang
e OUTPUT PARAMETERS
within [0..M-1, 0..Min(M,N)-1]. S - modified matrix
if UNeeded=2, U contains matrix U wholly. Array whose
indexes range within [0..M-1, 0..M-1]. NOTE 1: when S[i,j] is exactly zero after modification, it is deleted
VT - if VTNeeded=0, VT isn from the table.
are not calculated.
if VTNeeded=1, VT contains right singular vectors (firs -- ALGLIB PROJECT --
t Copyright 14.10.2011 by Bochkanov Sergey
min(M,N) rows of matrix V^T). Array whose indexes range *************************************************************************/
within [0..min(M,N)-1, 0..N-1]. void sparseadd(const sparsematrix &s, const ae_int_t i, const ae_int_t j, c
if VTNeeded=2, VT contains matrix V^T wholly. Array who onst double v);
se
indexes range within [0..N-1, 0..N-1]. /*************************************************************************
This function modifies S[i,j] - element of the sparse matrix. Matrix must
be in a Hash-Table mode.
In case new value of S[i,j] is zero, this element is deleted from the
table.
INPUT PARAMETERS
S - sparse M*N matrix in Hash-Table representation.
Exception will be thrown for CRS matrix.
I - row index of the element to modify, 0<=I<M
J - column index of the element to modify, 0<=J<N
V - value to set, must be finite number, can be zero
OUTPUT PARAMETERS
S - modified matrix
NOTE: this function has no effect when called with zero V for non-
existent element.
-- ALGLIB PROJECT --
Copyright 14.10.2011 by Bochkanov Sergey
*************************************************************************/
void sparseset(const sparsematrix &s, const ae_int_t i, const ae_int_t j, c
onst double v);
/*************************************************************************
This function returns S[i,j] - element of the sparse matrix. Matrix can
be in any mode (Hash-Table or CRS), but this function is less efficient
for CRS matrices. Hash-Table matrices can find element in O(1) time,
while CRS matrices need O(RS) time, where RS is an number of non-zero
elements in a row.
INPUT PARAMETERS
S - sparse M*N matrix in Hash-Table representation.
Exception will be thrown for CRS matrix.
I - row index of the element to modify, 0<=I<M
J - column index of the element to modify, 0<=J<N
RESULT
value of S[I,J] or zero (in case no element with such index is found)
-- ALGLIB PROJECT --
Copyright 14.10.2011 by Bochkanov Sergey
*************************************************************************/
double sparseget(const sparsematrix &s, const ae_int_t i, const ae_int_t j)
;
/*************************************************************************
This function converts matrix to CRS format.
Some algorithms (linear algebra ones, for example) require matrices in
CRS format.
INPUT PARAMETERS
S - sparse M*N matrix.
OUTPUT PARAMETERS
S - matrix in CRS format
NOTE: this function has no effect when called with matrix which is
already in CRS mode.
-- ALGLIB PROJECT --
Copyright 14.10.2011 by Bochkanov Sergey
*************************************************************************/
void sparseconverttocrs(const sparsematrix &s);
/*************************************************************************
This function calculates matrix-vector product S*x. Matrix S must be
stored in CRS format (exception will be thrown otherwise).
INPUT PARAMETERS
S - sparse M*N matrix in CRS format (you MUST convert it
to CRS before calling this function).
X - array[N], input vector. For performance reasons we
make only quick checks - we check that array size is
at least N, but we do not check for NAN's or INF's.
Y - output buffer, possibly preallocated. In case buffer
size is too small to store result, this buffer is
automatically resized.
OUTPUT PARAMETERS
Y - array[M], S*x
NOTE: this function throws exception when called for non-CRS matrix. You
must convert your matrix with SparseConvertToCRS() before using this
function.
-- ALGLIB PROJECT --
Copyright 14.10.2011 by Bochkanov Sergey
*************************************************************************/
void sparsemv(const sparsematrix &s, const real_1d_array &x, real_1d_array
&y);
/*************************************************************************
This function calculates matrix-vector product S^T*x. Matrix S must be
stored in CRS format (exception will be thrown otherwise).
INPUT PARAMETERS
S - sparse M*N matrix in CRS format (you MUST convert it
to CRS before calling this function).
X - array[M], input vector. For performance reasons we
make only quick checks - we check that array size is
at least M, but we do not check for NAN's or INF's.
Y - output buffer, possibly preallocated. In case buffer
size is too small to store result, this buffer is
automatically resized.
OUTPUT PARAMETERS
Y - array[N], S^T*x
NOTE: this function throws exception when called for non-CRS matrix. You
must convert your matrix with SparseConvertToCRS() before using this
function.
-- ALGLIB PROJECT --
Copyright 14.10.2011 by Bochkanov Sergey
*************************************************************************/
void sparsemtv(const sparsematrix &s, const real_1d_array &x, real_1d_array
&y);
/*************************************************************************
This function simultaneously calculates two matrix-vector products:
S*x and S^T*x.
S must be square (non-rectangular) matrix stored in CRS format (exception
will be thrown otherwise).
INPUT PARAMETERS
S - sparse N*N matrix in CRS format (you MUST convert it
to CRS before calling this function).
X - array[N], input vector. For performance reasons we
make only quick checks - we check that array size is
at least N, but we do not check for NAN's or INF's.
Y0 - output buffer, possibly preallocated. In case buffer
size is too small to store result, this buffer is
automatically resized.
Y1 - output buffer, possibly preallocated. In case buffer
size is too small to store result, this buffer is
automatically resized.
OUTPUT PARAMETERS
Y0 - array[N], S*x
Y1 - array[N], S^T*x
NOTE: this function throws exception when called for non-CRS matrix. You
must convert your matrix with SparseConvertToCRS() before using this
function. It also throws exception when S is non-square.
-- ALGLIB PROJECT --
Copyright 14.10.2011 by Bochkanov Sergey
*************************************************************************/
void sparsemv2(const sparsematrix &s, const real_1d_array &x, real_1d_array
&y0, real_1d_array &y1);
/*************************************************************************
This function calculates matrix-vector product S*x, when S is symmetric
matrix. Matrix S must be stored in CRS format (exception will be
thrown otherwise).
INPUT PARAMETERS
S - sparse M*M matrix in CRS format (you MUST convert it
to CRS before calling this function).
X - array[N], input vector. For performance reasons we
make only quick checks - we check that array size is
at least N, but we do not check for NAN's or INF's.
Y - output buffer, possibly preallocated. In case buffer
size is too small to store result, this buffer is
automatically resized.
OUTPUT PARAMETERS
Y - array[M], S*x
NOTE: this function throws exception when called for non-CRS matrix. You
must convert your matrix with SparseConvertToCRS() before using this
function.
-- ALGLIB PROJECT --
Copyright 14.10.2011 by Bochkanov Sergey
*************************************************************************/
void sparsesmv(const sparsematrix &s, const bool isupper, const real_1d_arr
ay &x, real_1d_array &y);
/*************************************************************************
This function calculates matrix-matrix product S*A. Matrix S must be
stored in CRS format (exception will be thrown otherwise).
INPUT PARAMETERS
S - sparse M*N matrix in CRS format (you MUST convert it
to CRS before calling this function).
A - array[N][K], input dense matrix. For performance reaso
ns
we make only quick checks - we check that array size
is at least N, but we do not check for NAN's or INF's.
K - number of columns of matrix (A).
B - output buffer, possibly preallocated. In case buffer
size is too small to store result, this buffer is
automatically resized.
OUTPUT PARAMETERS
B - array[M][K], S*A
NOTE: this function throws exception when called for non-CRS matrix. You
must convert your matrix with SparseConvertToCRS() before using this
function.
-- ALGLIB PROJECT --
Copyright 14.10.2011 by Bochkanov Sergey
*************************************************************************/
void sparsemm(const sparsematrix &s, const real_2d_array &a, const ae_int_t
k, real_2d_array &b);
/*************************************************************************
This function calculates matrix-matrix product S^T*A. Matrix S must be
stored in CRS format (exception will be thrown otherwise).
INPUT PARAMETERS
S - sparse M*N matrix in CRS format (you MUST convert it
to CRS before calling this function).
A - array[M][K], input dense matrix. For performance reason
s
we make only quick checks - we check that array size i
s
at least M, but we do not check for NAN's or INF's.
K - number of columns of matrix (A).
B - output buffer, possibly preallocated. In case buffer
size is too small to store result, this buffer is
automatically resized.
OUTPUT PARAMETERS
B - array[N][K], S^T*A
NOTE: this function throws exception when called for non-CRS matrix. You
must convert your matrix with SparseConvertToCRS() before using this
function.
-- ALGLIB PROJECT --
Copyright 14.10.2011 by Bochkanov Sergey
*************************************************************************/
void sparsemtm(const sparsematrix &s, const real_2d_array &a, const ae_int_
t k, real_2d_array &b);
/*************************************************************************
This function simultaneously calculates two matrix-matrix products:
S*A and S^T*A.
S must be square (non-rectangular) matrix stored in CRS format (exception
will be thrown otherwise).
INPUT PARAMETERS
S - sparse N*N matrix in CRS format (you MUST convert it
to CRS before calling this function).
A - array[N][K], input dense matrix. For performance reason
s
we make only quick checks - we check that array size i
s
at least N, but we do not check for NAN's or INF's.
K - number of columns of matrix (A).
B0 - output buffer, possibly preallocated. In case buffer
size is too small to store result, this buffer is
automatically resized.
B1 - output buffer, possibly preallocated. In case buffer
size is too small to store result, this buffer is
automatically resized.
OUTPUT PARAMETERS
B0 - array[N][K], S*A
B1 - array[N][K], S^T*A
NOTE: this function throws exception when called for non-CRS matrix. You
must convert your matrix with SparseConvertToCRS() before using this
function. It also throws exception when S is non-square.
-- ALGLIB PROJECT --
Copyright 14.10.2011 by Bochkanov Sergey
*************************************************************************/
void sparsemm2(const sparsematrix &s, const real_2d_array &a, const ae_int_
t k, real_2d_array &b0, real_2d_array &b1);
/*************************************************************************
This function calculates matrix-matrix product S*A, when S is symmetric
matrix. Matrix S must be stored in CRS format (exception will be
thrown otherwise).
INPUT PARAMETERS
S - sparse M*M matrix in CRS format (you MUST convert it
to CRS before calling this function).
A - array[N][K], input dense matrix. For performance reason
s
we make only quick checks - we check that array size is
at least N, but we do not check for NAN's or INF's.
K - number of columns of matrix (A).
B - output buffer, possibly preallocated. In case buffer
size is too small to store result, this buffer is
automatically resized.
OUTPUT PARAMETERS
B - array[M][K], S*A
NOTE: this function throws exception when called for non-CRS matrix. You
must convert your matrix with SparseConvertToCRS() before using this
function.
-- ALGLIB PROJECT --
Copyright 14.10.2011 by Bochkanov Sergey
*************************************************************************/
void sparsesmm(const sparsematrix &s, const bool isupper, const real_2d_arr
ay &a, const ae_int_t k, real_2d_array &b);
/*************************************************************************
This procedure resizes Hash-Table matrix. It can be called when you have
deleted too many elements from the matrix, and you want to free unneeded
memory.
-- ALGLIB PROJECT --
Copyright 14.10.2011 by Bochkanov Sergey
*************************************************************************/
void sparseresizematrix(const sparsematrix &s);
/*************************************************************************
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 -- -- ALGLIB --
Copyright 2005 by Bochkanov Sergey Copyright 06.12.2011 by Bochkanov Sergey
*************************************************************************/ *************************************************************************/
bool rmatrixsvd(const real_2d_array &a, const ae_int_t m, const ae_int_t n, void normestimatorcreate(const ae_int_t m, const ae_int_t n, const ae_int_t
const ae_int_t uneeded, const ae_int_t vtneeded, const ae_int_t additional nstart, const ae_int_t nits, normestimatorstate &state);
memory, real_1d_array &w, real_2d_array &u, real_2d_array &vt);
/*************************************************************************
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
*************************************************************************/
void normestimatorsetseed(const normestimatorstate &state, const ae_int_t s
eedval);
/*************************************************************************
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
*************************************************************************/
void normestimatorestimatesparse(const normestimatorstate &state, const spa
rsematrix &a);
/*************************************************************************
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
*************************************************************************/
void normestimatorresults(const normestimatorstate &state, double &nrm);
/************************************************************************* /*************************************************************************
Determinant calculation of the matrix given by its LU decomposition. Determinant calculation of the matrix given by its LU decomposition.
Input parameters: Input parameters:
A - LU decomposition of the matrix (output of A - LU decomposition of the matrix (output of
RMatrixLU subroutine). RMatrixLU subroutine).
Pivots - table of permutations which were made during Pivots - table of permutations which were made during
the LU decomposition. the LU decomposition.
Output of RMatrixLU subroutine. Output of RMatrixLU subroutine.
skipping to change at line 3465 skipping to change at line 4221
ae_int_t n, ae_int_t n,
/* Complex */ ae_vector* tau, /* Complex */ ae_vector* tau,
ae_int_t qrows, ae_int_t qrows,
/* Complex */ ae_matrix* q, /* Complex */ ae_matrix* q,
ae_state *_state); ae_state *_state);
void cmatrixlqunpackl(/* Complex */ ae_matrix* a, void cmatrixlqunpackl(/* Complex */ ae_matrix* a,
ae_int_t m, ae_int_t m,
ae_int_t n, ae_int_t n,
/* Complex */ ae_matrix* l, /* Complex */ ae_matrix* l,
ae_state *_state); ae_state *_state);
void rmatrixqrbasecase(/* Real */ ae_matrix* a,
ae_int_t m,
ae_int_t n,
/* Real */ ae_vector* work,
/* Real */ ae_vector* t,
/* Real */ ae_vector* tau,
ae_state *_state);
void rmatrixlqbasecase(/* Real */ ae_matrix* a,
ae_int_t m,
ae_int_t n,
/* Real */ ae_vector* work,
/* Real */ ae_vector* t,
/* Real */ ae_vector* tau,
ae_state *_state);
void rmatrixbd(/* Real */ ae_matrix* a, void rmatrixbd(/* Real */ ae_matrix* a,
ae_int_t m, ae_int_t m,
ae_int_t n, ae_int_t n,
/* Real */ ae_vector* tauq, /* Real */ ae_vector* tauq,
/* Real */ ae_vector* taup, /* Real */ ae_vector* taup,
ae_state *_state); ae_state *_state);
void rmatrixbdunpackq(/* Real */ ae_matrix* qp, void rmatrixbdunpackq(/* Real */ ae_matrix* qp,
ae_int_t m, ae_int_t m,
ae_int_t n, ae_int_t n,
/* Real */ ae_vector* tauq, /* Real */ ae_vector* tauq,
skipping to change at line 3551 skipping to change at line 4321
/* Complex */ ae_vector* tau, /* Complex */ ae_vector* tau,
/* Real */ ae_vector* d, /* Real */ ae_vector* d,
/* Real */ ae_vector* e, /* Real */ ae_vector* e,
ae_state *_state); ae_state *_state);
void hmatrixtdunpackq(/* Complex */ ae_matrix* a, void hmatrixtdunpackq(/* Complex */ ae_matrix* a,
ae_int_t n, ae_int_t n,
ae_bool isupper, ae_bool isupper,
/* Complex */ ae_vector* tau, /* Complex */ ae_vector* tau,
/* Complex */ ae_matrix* q, /* Complex */ ae_matrix* q,
ae_state *_state); ae_state *_state);
ae_bool rmatrixbdsvd(/* Real */ ae_vector* d,
/* Real */ ae_vector* e,
ae_int_t n,
ae_bool isupper,
ae_bool isfractionalaccuracyrequired,
/* Real */ ae_matrix* u,
ae_int_t nru,
/* Real */ ae_matrix* c,
ae_int_t ncc,
/* Real */ ae_matrix* vt,
ae_int_t ncvt,
ae_state *_state);
ae_bool bidiagonalsvddecomposition(/* Real */ ae_vector* d,
/* Real */ ae_vector* e,
ae_int_t n,
ae_bool isupper,
ae_bool isfractionalaccuracyrequired,
/* Real */ ae_matrix* u,
ae_int_t nru,
/* Real */ ae_matrix* c,
ae_int_t ncc,
/* Real */ ae_matrix* vt,
ae_int_t ncvt,
ae_state *_state);
ae_bool rmatrixsvd(/* Real */ ae_matrix* a,
ae_int_t m,
ae_int_t n,
ae_int_t uneeded,
ae_int_t vtneeded,
ae_int_t additionalmemory,
/* Real */ ae_vector* w,
/* Real */ ae_matrix* u,
/* Real */ ae_matrix* vt,
ae_state *_state);
ae_bool smatrixevd(/* Real */ ae_matrix* a, ae_bool smatrixevd(/* Real */ ae_matrix* a,
ae_int_t n, ae_int_t n,
ae_int_t zneeded, ae_int_t zneeded,
ae_bool isupper, ae_bool isupper,
/* Real */ ae_vector* d, /* Real */ ae_vector* d,
/* Real */ ae_matrix* z, /* Real */ ae_matrix* z,
ae_state *_state); ae_state *_state);
ae_bool smatrixevdr(/* Real */ ae_matrix* a, ae_bool smatrixevdr(/* Real */ ae_matrix* a,
ae_int_t n, ae_int_t n,
ae_int_t zneeded, ae_int_t zneeded,
skipping to change at line 3854 skipping to change at line 4658
void cmatrixtrinverse(/* Complex */ ae_matrix* a, void cmatrixtrinverse(/* Complex */ ae_matrix* a,
ae_int_t n, ae_int_t n,
ae_bool isupper, ae_bool isupper,
ae_bool isunit, ae_bool isunit,
ae_int_t* info, ae_int_t* info,
matinvreport* rep, matinvreport* rep,
ae_state *_state); ae_state *_state);
ae_bool _matinvreport_init(matinvreport* p, ae_state *_state, ae_bool make_ automatic); ae_bool _matinvreport_init(matinvreport* p, ae_state *_state, ae_bool make_ automatic);
ae_bool _matinvreport_init_copy(matinvreport* dst, matinvreport* src, ae_st ate *_state, ae_bool make_automatic); ae_bool _matinvreport_init_copy(matinvreport* dst, matinvreport* src, ae_st ate *_state, ae_bool make_automatic);
void _matinvreport_clear(matinvreport* p); void _matinvreport_clear(matinvreport* p);
ae_bool rmatrixbdsvd(/* Real */ ae_vector* d,
/* Real */ ae_vector* e,
ae_int_t n,
ae_bool isupper,
ae_bool isfractionalaccuracyrequired,
/* Real */ ae_matrix* u,
ae_int_t nru,
/* Real */ ae_matrix* c,
ae_int_t ncc,
/* Real */ ae_matrix* vt,
ae_int_t ncvt,
ae_state *_state);
ae_bool bidiagonalsvddecomposition(/* Real */ ae_vector* d,
/* Real */ ae_vector* e,
ae_int_t n,
ae_bool isupper,
ae_bool isfractionalaccuracyrequired,
/* Real */ ae_matrix* u,
ae_int_t nru,
/* Real */ ae_matrix* c,
ae_int_t ncc,
/* Real */ ae_matrix* vt,
ae_int_t ncvt,
ae_state *_state);
ae_bool rmatrixsvd(/* Real */ ae_matrix* a,
ae_int_t m,
ae_int_t n,
ae_int_t uneeded,
ae_int_t vtneeded,
ae_int_t additionalmemory,
/* Real */ ae_vector* w,
/* Real */ ae_matrix* u,
/* Real */ ae_matrix* vt,
ae_state *_state);
void fblscholeskysolve(/* Real */ ae_matrix* cha, void fblscholeskysolve(/* Real */ ae_matrix* cha,
double sqrtscalea, double sqrtscalea,
ae_int_t n, ae_int_t n,
ae_bool isupper, ae_bool isupper,
/* Real */ ae_vector* xb, /* Real */ ae_vector* xb,
/* Real */ ae_vector* tmp, /* Real */ ae_vector* tmp,
ae_state *_state); ae_state *_state);
void fblssolvecgx(/* Real */ ae_matrix* a, void fblssolvecgx(/* Real */ ae_matrix* a,
ae_int_t m, ae_int_t m,
ae_int_t n, ae_int_t n,
skipping to change at line 3909 skipping to change at line 4679
/* Real */ ae_vector* b, /* Real */ ae_vector* b,
/* Real */ ae_vector* x, /* Real */ ae_vector* x,
/* Real */ ae_vector* buf, /* Real */ ae_vector* buf,
ae_state *_state); ae_state *_state);
void fblscgcreate(/* Real */ ae_vector* x, void fblscgcreate(/* Real */ ae_vector* x,
/* Real */ ae_vector* b, /* Real */ ae_vector* b,
ae_int_t n, ae_int_t n,
fblslincgstate* state, fblslincgstate* state,
ae_state *_state); ae_state *_state);
ae_bool fblscgiteration(fblslincgstate* state, ae_state *_state); ae_bool fblscgiteration(fblslincgstate* state, ae_state *_state);
void fblssolvels(/* Real */ ae_matrix* a,
/* Real */ ae_vector* b,
ae_int_t m,
ae_int_t n,
/* Real */ ae_vector* tmp0,
/* Real */ ae_vector* tmp1,
/* Real */ ae_vector* tmp2,
ae_state *_state);
ae_bool _fblslincgstate_init(fblslincgstate* p, ae_state *_state, ae_bool m ake_automatic); ae_bool _fblslincgstate_init(fblslincgstate* p, ae_state *_state, ae_bool m ake_automatic);
ae_bool _fblslincgstate_init_copy(fblslincgstate* dst, fblslincgstate* src, ae_state *_state, ae_bool make_automatic); ae_bool _fblslincgstate_init_copy(fblslincgstate* dst, fblslincgstate* src, ae_state *_state, ae_bool make_automatic);
void _fblslincgstate_clear(fblslincgstate* p); void _fblslincgstate_clear(fblslincgstate* p);
void sparsecreate(ae_int_t m,
ae_int_t n,
ae_int_t k,
sparsematrix* s,
ae_state *_state);
void sparsecreatecrs(ae_int_t m,
ae_int_t n,
/* Integer */ ae_vector* ner,
sparsematrix* s,
ae_state *_state);
void sparsecopy(sparsematrix* s0, sparsematrix* s1, ae_state *_state);
void sparseadd(sparsematrix* s,
ae_int_t i,
ae_int_t j,
double v,
ae_state *_state);
void sparseset(sparsematrix* s,
ae_int_t i,
ae_int_t j,
double v,
ae_state *_state);
double sparseget(sparsematrix* s,
ae_int_t i,
ae_int_t j,
ae_state *_state);
void sparseconverttocrs(sparsematrix* s, ae_state *_state);
void sparsemv(sparsematrix* s,
/* Real */ ae_vector* x,
/* Real */ ae_vector* y,
ae_state *_state);
void sparsemtv(sparsematrix* s,
/* Real */ ae_vector* x,
/* Real */ ae_vector* y,
ae_state *_state);
void sparsemv2(sparsematrix* s,
/* Real */ ae_vector* x,
/* Real */ ae_vector* y0,
/* Real */ ae_vector* y1,
ae_state *_state);
void sparsesmv(sparsematrix* s,
ae_bool isupper,
/* Real */ ae_vector* x,
/* Real */ ae_vector* y,
ae_state *_state);
void sparsemm(sparsematrix* s,
/* Real */ ae_matrix* a,
ae_int_t k,
/* Real */ ae_matrix* b,
ae_state *_state);
void sparsemtm(sparsematrix* s,
/* Real */ ae_matrix* a,
ae_int_t k,
/* Real */ ae_matrix* b,
ae_state *_state);
void sparsemm2(sparsematrix* s,
/* Real */ ae_matrix* a,
ae_int_t k,
/* Real */ ae_matrix* b0,
/* Real */ ae_matrix* b1,
ae_state *_state);
void sparsesmm(sparsematrix* s,
ae_bool isupper,
/* Real */ ae_matrix* a,
ae_int_t k,
/* Real */ ae_matrix* b,
ae_state *_state);
void sparseresizematrix(sparsematrix* s, ae_state *_state);
double sparsegetaveragelengthofchain(sparsematrix* s, ae_state *_state);
ae_bool _sparsematrix_init(sparsematrix* p, ae_state *_state, ae_bool make_
automatic);
ae_bool _sparsematrix_init_copy(sparsematrix* dst, sparsematrix* src, ae_st
ate *_state, ae_bool make_automatic);
void _sparsematrix_clear(sparsematrix* p);
void normestimatorcreate(ae_int_t m,
ae_int_t n,
ae_int_t nstart,
ae_int_t nits,
normestimatorstate* state,
ae_state *_state);
void normestimatorsetseed(normestimatorstate* state,
ae_int_t seedval,
ae_state *_state);
ae_bool normestimatoriteration(normestimatorstate* state,
ae_state *_state);
void normestimatorestimatesparse(normestimatorstate* state,
sparsematrix* a,
ae_state *_state);
void normestimatorresults(normestimatorstate* state,
double* nrm,
ae_state *_state);
void normestimatorrestart(normestimatorstate* state, ae_state *_state);
ae_bool _normestimatorstate_init(normestimatorstate* p, ae_state *_state, a
e_bool make_automatic);
ae_bool _normestimatorstate_init_copy(normestimatorstate* dst, normestimato
rstate* src, ae_state *_state, ae_bool make_automatic);
void _normestimatorstate_clear(normestimatorstate* p);
double rmatrixludet(/* Real */ ae_matrix* a, double rmatrixludet(/* Real */ ae_matrix* a,
/* Integer */ ae_vector* pivots, /* Integer */ ae_vector* pivots,
ae_int_t n, ae_int_t n,
ae_state *_state); ae_state *_state);
double rmatrixdet(/* Real */ ae_matrix* a, double rmatrixdet(/* Real */ ae_matrix* a,
ae_int_t n, ae_int_t n,
ae_state *_state); ae_state *_state);
ae_complex cmatrixludet(/* Complex */ ae_matrix* a, ae_complex cmatrixludet(/* Complex */ ae_matrix* a,
/* Integer */ ae_vector* pivots, /* Integer */ ae_vector* pivots,
ae_int_t n, ae_int_t n,
 End of changes. 43 change blocks. 
169 lines changed or deleted 1068 lines changed or added


 optimization.h   optimization.h 
skipping to change at line 23 skipping to change at line 23
GNU General Public License for more details. GNU General Public License for more details.
A copy of the GNU General Public License is available at A copy of the GNU General Public License is available at
http://www.fsf.org/licensing/licenses http://www.fsf.org/licensing/licenses
>>> END OF LICENSE >>> >>> END OF LICENSE >>>
*************************************************************************/ *************************************************************************/
#ifndef _optimization_pkg_h #ifndef _optimization_pkg_h
#define _optimization_pkg_h #define _optimization_pkg_h
#include "ap.h" #include "ap.h"
#include "alglibinternal.h" #include "alglibinternal.h"
#include "alglibmisc.h"
#include "linalg.h" #include "linalg.h"
#include "alglibmisc.h"
#include "solvers.h" #include "solvers.h"
///////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////
// //
// THIS SECTION CONTAINS COMPUTATIONAL CORE DECLARATIONS (DATATYPES) // THIS SECTION CONTAINS COMPUTATIONAL CORE DECLARATIONS (DATATYPES)
// //
///////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////
namespace alglib_impl namespace alglib_impl
{ {
typedef struct typedef struct
skipping to change at line 65 skipping to change at line 65
ae_int_t k; ae_int_t k;
ae_vector xk; ae_vector xk;
ae_vector dk; ae_vector dk;
ae_vector xn; ae_vector xn;
ae_vector dn; ae_vector dn;
ae_vector d; ae_vector d;
double fold; double fold;
double stp; double stp;
double curstpmax; double curstpmax;
ae_vector yk; ae_vector yk;
double laststep; double lastgoodstep;
double lastscaledstep; double lastscaledstep;
ae_int_t mcinfo; ae_int_t mcinfo;
ae_bool innerresetneeded; ae_bool innerresetneeded;
ae_bool terminationneeded; ae_bool terminationneeded;
double trimthreshold; double trimthreshold;
ae_int_t rstimer; ae_int_t rstimer;
ae_vector x; ae_vector x;
double f; double f;
ae_vector g; ae_vector g;
ae_bool needf; ae_bool needf;
skipping to change at line 134 skipping to change at line 134
ae_bool xupdated; ae_bool xupdated;
rcommstate rstate; rcommstate rstate;
ae_int_t repinneriterationscount; ae_int_t repinneriterationscount;
ae_int_t repouteriterationscount; ae_int_t repouteriterationscount;
ae_int_t repnfev; ae_int_t repnfev;
ae_int_t repterminationtype; ae_int_t repterminationtype;
double repdebugeqerr; double repdebugeqerr;
double repdebugfs; double repdebugfs;
double repdebugff; double repdebugff;
double repdebugdx; double repdebugdx;
ae_int_t repdebugfeasqpits;
ae_int_t repdebugfeasgpaits;
ae_vector xcur; ae_vector xcur;
ae_vector xprev; ae_vector xprev;
ae_vector xstart; ae_vector xstart;
ae_int_t itsleft; ae_int_t itsleft;
ae_vector xend; ae_vector xend;
ae_vector lastg; ae_vector lastg;
double trimthreshold; double trimthreshold;
ae_matrix ceoriginal; ae_matrix ceoriginal;
ae_matrix ceeffective; ae_matrix ceeffective;
ae_matrix cecurrent; ae_matrix cecurrent;
skipping to change at line 196 skipping to change at line 198
typedef struct typedef struct
{ {
ae_int_t inneriterationscount; ae_int_t inneriterationscount;
ae_int_t outeriterationscount; ae_int_t outeriterationscount;
ae_int_t nfev; ae_int_t nfev;
ae_int_t terminationtype; ae_int_t terminationtype;
double debugeqerr; double debugeqerr;
double debugfs; double debugfs;
double debugff; double debugff;
double debugdx; double debugdx;
ae_int_t debugfeasqpits;
ae_int_t debugfeasgpaits;
} minbleicreport; } minbleicreport;
typedef struct typedef struct
{ {
ae_int_t n; ae_int_t n;
ae_int_t m; ae_int_t m;
double epsg; double epsg;
double epsf; double epsf;
double epsx; double epsx;
ae_int_t maxits; ae_int_t maxits;
ae_bool xrep; ae_bool xrep;
skipping to change at line 257 skipping to change at line 261
ae_int_t iterationscount; ae_int_t iterationscount;
ae_int_t nfev; ae_int_t nfev;
ae_int_t terminationtype; ae_int_t terminationtype;
} minlbfgsreport; } minlbfgsreport;
typedef struct typedef struct
{ {
ae_int_t n; ae_int_t n;
ae_int_t algokind; ae_int_t algokind;
ae_int_t akind; ae_int_t akind;
ae_matrix densea; ae_matrix densea;
ae_vector diaga;
ae_vector b; ae_vector b;
ae_vector bndl; ae_vector bndl;
ae_vector bndu; ae_vector bndu;
ae_vector havebndl; ae_vector havebndl;
ae_vector havebndu; ae_vector havebndu;
ae_vector xorigin; ae_vector xorigin;
ae_vector startx; ae_vector startx;
ae_bool havex; ae_bool havex;
ae_vector xc; ae_vector xc;
ae_vector gc; ae_vector gc;
ae_vector activeconstraints; ae_vector dc;
ae_vector prevactiveconstraints; ae_vector xcand0;
ae_vector xcand1;
ae_vector xn;
ae_vector pg;
ae_vector activeset;
double constterm; double constterm;
ae_vector workbndl; ae_vector workbndl;
ae_vector workbndu; ae_vector workbndu;
ae_int_t repinneriterationscount; ae_int_t repinneriterationscount;
ae_int_t repouteriterationscount; ae_int_t repouteriterationscount;
ae_int_t repncholesky; ae_int_t repncholesky;
ae_int_t repnmv; ae_int_t repnmv;
ae_int_t repterminationtype; ae_int_t repterminationtype;
ae_vector tmp0; ae_vector tmp0;
ae_vector tmp1; ae_vector tmp1;
skipping to change at line 595 skipping to change at line 602
minbleicreport& operator=(const minbleicreport &rhs); minbleicreport& operator=(const minbleicreport &rhs);
virtual ~minbleicreport(); virtual ~minbleicreport();
ae_int_t &inneriterationscount; ae_int_t &inneriterationscount;
ae_int_t &outeriterationscount; ae_int_t &outeriterationscount;
ae_int_t &nfev; ae_int_t &nfev;
ae_int_t &terminationtype; ae_int_t &terminationtype;
double &debugeqerr; double &debugeqerr;
double &debugfs; double &debugfs;
double &debugff; double &debugff;
double &debugdx; double &debugdx;
ae_int_t &debugfeasqpits;
ae_int_t &debugfeasgpaits;
}; };
/************************************************************************* /*************************************************************************
*************************************************************************/ *************************************************************************/
class _minlbfgsstate_owner class _minlbfgsstate_owner
{ {
public: public:
_minlbfgsstate_owner(); _minlbfgsstate_owner();
skipping to change at line 701 skipping to change at line 710
process are counted) process are counted)
* TerminationType completion code (see below) * TerminationType completion code (see below)
Completion codes: Completion codes:
* -5 inappropriate solver was used: * -5 inappropriate solver was used:
* Cholesky solver for semidefinite or indefinite problems * Cholesky solver for semidefinite or indefinite problems
* Cholesky solver for problems with non-boundary constraints * Cholesky solver for problems with non-boundary constraints
* -3 inconsistent constraints (or, maybe, feasible point is * -3 inconsistent constraints (or, maybe, feasible point is
too hard to find). If you are sure that constraints are feasible, too hard to find). If you are sure that constraints are feasible,
try to restart optimizer with better initial approximation. try to restart optimizer with better initial approximation.
* -1 solver error
* 4 successful completion * 4 successful completion
* 5 MaxIts steps was taken * 5 MaxIts steps was taken
* 7 stopping conditions are too stringent, * 7 stopping conditions are too stringent,
further improvement is impossible, further improvement is impossible,
X contains best point found so far. X contains best point found so far.
*************************************************************************/ *************************************************************************/
class _minqpreport_owner class _minqpreport_owner
{ {
public: public:
_minqpreport_owner(); _minqpreport_owner();
skipping to change at line 3068 skipping to change at line 3078
void minasarestartfrom(const minasastate &state, const real_1d_array &x, co nst real_1d_array &bndl, const real_1d_array &bndu); void minasarestartfrom(const minasastate &state, const real_1d_array &x, co nst real_1d_array &bndl, const real_1d_array &bndu);
} }
///////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////
// //
// THIS SECTION CONTAINS COMPUTATIONAL CORE DECLARATIONS (FUNCTIONS) // THIS SECTION CONTAINS COMPUTATIONAL CORE DECLARATIONS (FUNCTIONS)
// //
///////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////
namespace alglib_impl namespace alglib_impl
{ {
void trimprepare(double f, double* threshold, ae_state *_state);
void trimfunction(double* f,
/* Real */ ae_vector* g,
ae_int_t n,
double threshold,
ae_state *_state);
ae_bool enforceboundaryconstraints(/* Real */ ae_vector* x,
/* Real */ ae_vector* bl,
/* Boolean */ ae_vector* havebl,
/* Real */ ae_vector* bu,
/* Boolean */ ae_vector* havebu,
ae_int_t nmain,
ae_int_t nslack,
ae_state *_state);
void projectgradientintobc(/* Real */ ae_vector* x,
/* Real */ ae_vector* g,
/* Real */ ae_vector* bl,
/* Boolean */ ae_vector* havebl,
/* Real */ ae_vector* bu,
/* Boolean */ ae_vector* havebu,
ae_int_t nmain,
ae_int_t nslack,
ae_state *_state);
void calculatestepbound(/* Real */ ae_vector* x,
/* Real */ ae_vector* d,
double alpha,
/* Real */ ae_vector* bndl,
/* Boolean */ ae_vector* havebndl,
/* Real */ ae_vector* bndu,
/* Boolean */ ae_vector* havebndu,
ae_int_t nmain,
ae_int_t nslack,
ae_int_t* variabletofreeze,
double* valuetofreeze,
double* maxsteplen,
ae_state *_state);
ae_int_t postprocessboundedstep(/* Real */ ae_vector* x,
/* Real */ ae_vector* xprev,
/* Real */ ae_vector* bndl,
/* Boolean */ ae_vector* havebndl,
/* Real */ ae_vector* bndu,
/* Boolean */ ae_vector* havebndu,
ae_int_t nmain,
ae_int_t nslack,
ae_int_t variabletofreeze,
double valuetofreeze,
double steptaken,
double maxsteplen,
ae_state *_state);
void filterdirection(/* Real */ ae_vector* d,
/* Real */ ae_vector* x,
/* Real */ ae_vector* bndl,
/* Boolean */ ae_vector* havebndl,
/* Real */ ae_vector* bndu,
/* Boolean */ ae_vector* havebndu,
/* Real */ ae_vector* s,
ae_int_t nmain,
ae_int_t nslack,
double droptol,
ae_state *_state);
ae_int_t numberofchangedconstraints(/* Real */ ae_vector* x,
/* Real */ ae_vector* xprev,
/* Real */ ae_vector* bndl,
/* Boolean */ ae_vector* havebndl,
/* Real */ ae_vector* bndu,
/* Boolean */ ae_vector* havebndu,
ae_int_t nmain,
ae_int_t nslack,
ae_state *_state);
ae_bool findfeasiblepoint(/* Real */ ae_vector* x,
/* Real */ ae_vector* bndl,
/* Boolean */ ae_vector* havebndl,
/* Real */ ae_vector* bndu,
/* Boolean */ ae_vector* havebndu,
ae_int_t nmain,
ae_int_t nslack,
/* Real */ ae_matrix* ce,
ae_int_t k,
double epsi,
ae_int_t* qpits,
ae_int_t* gpaits,
ae_state *_state);
void mincgcreate(ae_int_t n, void mincgcreate(ae_int_t n,
/* Real */ ae_vector* x, /* Real */ ae_vector* x,
mincgstate* state, mincgstate* state,
ae_state *_state); ae_state *_state);
void mincgcreatef(ae_int_t n, void mincgcreatef(ae_int_t n,
/* Real */ ae_vector* x, /* Real */ ae_vector* x,
double diffstep, double diffstep,
mincgstate* state, mincgstate* state,
ae_state *_state); ae_state *_state);
void mincgsetcond(mincgstate* state, void mincgsetcond(mincgstate* state,
 End of changes. 10 change blocks. 
5 lines changed or deleted 97 lines changed or added


 solvers.h   solvers.h 
skipping to change at line 47 skipping to change at line 47
} densesolverreport; } densesolverreport;
typedef struct typedef struct
{ {
double r2; double r2;
ae_matrix cx; ae_matrix cx;
ae_int_t n; ae_int_t n;
ae_int_t k; ae_int_t k;
} densesolverlsreport; } densesolverlsreport;
typedef struct typedef struct
{ {
normestimatorstate nes;
ae_vector rx;
ae_vector b;
ae_int_t n;
ae_int_t m;
ae_vector ui;
ae_vector uip1;
ae_vector vi;
ae_vector vip1;
ae_vector omegai;
ae_vector omegaip1;
double alphai;
double alphaip1;
double betai;
double betaip1;
double phibari;
double phibarip1;
double phii;
double rhobari;
double rhobarip1;
double rhoi;
double ci;
double si;
double theta;
double lambdai;
ae_vector d;
double anorm;
double bnorm2;
double dnorm;
double r2;
ae_vector x;
ae_vector mv;
ae_vector mtv;
double epsa;
double epsb;
double epsc;
ae_int_t maxits;
ae_bool xrep;
ae_bool xupdated;
ae_bool needmv;
ae_bool needmtv;
ae_bool needmv2;
ae_bool needvmv;
ae_bool needprec;
ae_int_t repiterationscount;
ae_int_t repnmv;
ae_int_t repterminationtype;
ae_bool running;
rcommstate rstate;
} linlsqrstate;
typedef struct
{
ae_int_t iterationscount;
ae_int_t nmv;
ae_int_t terminationtype;
} linlsqrreport;
typedef struct
{
ae_vector rx;
ae_vector b;
ae_int_t n;
ae_vector cx;
ae_vector cr;
ae_vector cz;
ae_vector p;
ae_vector r;
ae_vector z;
double alpha;
double beta;
double r2;
double meritfunction;
ae_vector x;
ae_vector mv;
ae_vector pv;
double vmv;
ae_vector startx;
double epsf;
ae_int_t maxits;
ae_int_t itsbeforerestart;
ae_int_t itsbeforerupdate;
ae_bool xrep;
ae_bool xupdated;
ae_bool needmv;
ae_bool needmtv;
ae_bool needmv2;
ae_bool needvmv;
ae_bool needprec;
ae_int_t repiterationscount;
ae_int_t repnmv;
ae_int_t repterminationtype;
ae_bool running;
rcommstate rstate;
} lincgstate;
typedef struct
{
ae_int_t iterationscount;
ae_int_t nmv;
ae_int_t terminationtype;
double r2;
} lincgreport;
typedef struct
{
ae_int_t n; ae_int_t n;
ae_int_t m; ae_int_t m;
double epsf; double epsf;
ae_int_t maxits; ae_int_t maxits;
ae_bool xrep; ae_bool xrep;
double stpmax; double stpmax;
ae_vector x; ae_vector x;
double f; double f;
ae_vector fi; ae_vector fi;
ae_matrix j; ae_matrix j;
skipping to change at line 147 skipping to change at line 249
densesolverlsreport& operator=(const densesolverlsreport &rhs); densesolverlsreport& operator=(const densesolverlsreport &rhs);
virtual ~densesolverlsreport(); virtual ~densesolverlsreport();
double &r2; double &r2;
real_2d_array cx; real_2d_array cx;
ae_int_t &n; ae_int_t &n;
ae_int_t &k; ae_int_t &k;
}; };
/************************************************************************* /*************************************************************************
This object stores state of the LinLSQR method.
You should use ALGLIB functions to work with this object.
*************************************************************************/
class _linlsqrstate_owner
{
public:
_linlsqrstate_owner();
_linlsqrstate_owner(const _linlsqrstate_owner &rhs);
_linlsqrstate_owner& operator=(const _linlsqrstate_owner &rhs);
virtual ~_linlsqrstate_owner();
alglib_impl::linlsqrstate* c_ptr();
alglib_impl::linlsqrstate* c_ptr() const;
protected:
alglib_impl::linlsqrstate *p_struct;
};
class linlsqrstate : public _linlsqrstate_owner
{
public:
linlsqrstate();
linlsqrstate(const linlsqrstate &rhs);
linlsqrstate& operator=(const linlsqrstate &rhs);
virtual ~linlsqrstate();
};
/*************************************************************************
*************************************************************************/
class _linlsqrreport_owner
{
public:
_linlsqrreport_owner();
_linlsqrreport_owner(const _linlsqrreport_owner &rhs);
_linlsqrreport_owner& operator=(const _linlsqrreport_owner &rhs);
virtual ~_linlsqrreport_owner();
alglib_impl::linlsqrreport* c_ptr();
alglib_impl::linlsqrreport* c_ptr() const;
protected:
alglib_impl::linlsqrreport *p_struct;
};
class linlsqrreport : public _linlsqrreport_owner
{
public:
linlsqrreport();
linlsqrreport(const linlsqrreport &rhs);
linlsqrreport& operator=(const linlsqrreport &rhs);
virtual ~linlsqrreport();
ae_int_t &iterationscount;
ae_int_t &nmv;
ae_int_t &terminationtype;
};
/*************************************************************************
This object stores state of the linear CG method.
You should use ALGLIB functions to work with this object.
Never try to access its fields directly!
*************************************************************************/
class _lincgstate_owner
{
public:
_lincgstate_owner();
_lincgstate_owner(const _lincgstate_owner &rhs);
_lincgstate_owner& operator=(const _lincgstate_owner &rhs);
virtual ~_lincgstate_owner();
alglib_impl::lincgstate* c_ptr();
alglib_impl::lincgstate* c_ptr() const;
protected:
alglib_impl::lincgstate *p_struct;
};
class lincgstate : public _lincgstate_owner
{
public:
lincgstate();
lincgstate(const lincgstate &rhs);
lincgstate& operator=(const lincgstate &rhs);
virtual ~lincgstate();
};
/*************************************************************************
*************************************************************************/
class _lincgreport_owner
{
public:
_lincgreport_owner();
_lincgreport_owner(const _lincgreport_owner &rhs);
_lincgreport_owner& operator=(const _lincgreport_owner &rhs);
virtual ~_lincgreport_owner();
alglib_impl::lincgreport* c_ptr();
alglib_impl::lincgreport* c_ptr() const;
protected:
alglib_impl::lincgreport *p_struct;
};
class lincgreport : public _lincgreport_owner
{
public:
lincgreport();
lincgreport(const lincgreport &rhs);
lincgreport& operator=(const lincgreport &rhs);
virtual ~lincgreport();
ae_int_t &iterationscount;
ae_int_t &nmv;
ae_int_t &terminationtype;
double &r2;
};
/*************************************************************************
*************************************************************************/ *************************************************************************/
class _nleqstate_owner class _nleqstate_owner
{ {
public: public:
_nleqstate_owner(); _nleqstate_owner();
_nleqstate_owner(const _nleqstate_owner &rhs); _nleqstate_owner(const _nleqstate_owner &rhs);
_nleqstate_owner& operator=(const _nleqstate_owner &rhs); _nleqstate_owner& operator=(const _nleqstate_owner &rhs);
virtual ~_nleqstate_owner(); virtual ~_nleqstate_owner();
alglib_impl::nleqstate* c_ptr(); alglib_impl::nleqstate* c_ptr();
skipping to change at line 878 skipping to change at line 1092
* K dim(Null(A)) * K dim(Null(A))
* CX array[0..N-1,0..K-1], kernel of A. * CX array[0..N-1,0..K-1], kernel of A.
Columns of CX store such vectors that A*CX[i]=0. Columns of CX store such vectors that A*CX[i]=0.
-- ALGLIB -- -- ALGLIB --
Copyright 24.08.2009 by Bochkanov Sergey Copyright 24.08.2009 by Bochkanov Sergey
*************************************************************************/ *************************************************************************/
void rmatrixsolvels(const real_2d_array &a, const ae_int_t nrows, const ae_ int_t ncols, const real_1d_array &b, const double threshold, ae_int_t &info , densesolverlsreport &rep, real_1d_array &x); void rmatrixsolvels(const real_2d_array &a, const ae_int_t nrows, const ae_ int_t ncols, const real_1d_array &b, const double threshold, ae_int_t &info , densesolverlsreport &rep, real_1d_array &x);
/************************************************************************* /*************************************************************************
This function initializes linear LSQR Solver. This solver is used to solve
non-symmetric (and, possibly, non-square) problems. Least squares solution
is returned for non-compatible systems.
USAGE:
1. User initializes algorithm state with LinLSQRCreate() call
2. User tunes solver parameters with LinLSQRSetCond() and other functions
3. User calls LinLSQRSolveSparse() function which takes algorithm state
and SparseMatrix object.
4. User calls LinLSQRResults() to get solution
5. Optionally, user may call LinLSQRSolveSparse() again to solve another
problem with different matrix and/or right part without reinitializing
LinLSQRState structure.
INPUT PARAMETERS:
M - number of rows in A
N - number of variables, N>0
OUTPUT PARAMETERS:
State - structure which stores algorithm state
-- ALGLIB --
Copyright 30.11.2011 by Bochkanov Sergey
*************************************************************************/
void linlsqrcreate(const ae_int_t m, const ae_int_t n, linlsqrstate &state)
;
/*************************************************************************
This function sets optional Tikhonov regularization coefficient.
It is zero by default.
INPUT PARAMETERS:
LambdaI - regularization factor, LambdaI>=0
OUTPUT PARAMETERS:
State - structure which stores algorithm state
-- ALGLIB --
Copyright 30.11.2011 by Bochkanov Sergey
*************************************************************************/
void linlsqrsetlambdai(const linlsqrstate &state, const double lambdai);
/*************************************************************************
Procedure for solution of A*x=b with sparse A.
INPUT PARAMETERS:
State - algorithm state
A - sparse M*N matrix in the CRS format (you MUST contvert it
to CRS format by calling SparseConvertToCRS() function
BEFORE you pass it to this function).
B - right part, array[M]
RESULT:
This function returns no result.
You can get solution by calling LinCGResults()
-- ALGLIB --
Copyright 30.11.2011 by Bochkanov Sergey
*************************************************************************/
void linlsqrsolvesparse(const linlsqrstate &state, const sparsematrix &a, c
onst real_1d_array &b);
/*************************************************************************
This function sets stopping criteria.
INPUT PARAMETERS:
EpsA - algorithm will be stopped if ||A^T*Rk||/(||A||*||Rk||)<=Eps
A.
EpsB - algorithm will be stopped if ||Rk||<=EpsB*||B||
MaxIts - algorithm will be stopped if number of iterations
more than MaxIts.
OUTPUT PARAMETERS:
State - structure which stores algorithm state
NOTE: if EpsA,EpsB,EpsC and MaxIts are zero then these variables will
be setted as default values.
-- ALGLIB --
Copyright 30.11.2011 by Bochkanov Sergey
*************************************************************************/
void linlsqrsetcond(const linlsqrstate &state, const double epsa, const dou
ble epsb, const ae_int_t maxits);
/*************************************************************************
LSQR solver: results.
This function must be called after LinLSQRSolve
INPUT PARAMETERS:
State - algorithm state
OUTPUT PARAMETERS:
X - array[N], solution
Rep - optimization report:
* Rep.TerminationType completetion code:
* 1 ||Rk||<=EpsB*||B||
* 4 ||A^T*Rk||/(||A||*||Rk||)<=EpsA
* 5 MaxIts steps was taken
* 7 rounding errors prevent further progress,
X contains best point found so far.
(sometimes returned on singular systems)
* Rep.IterationsCount contains iterations count
* NMV countains number of matrix-vector calculations
-- ALGLIB --
Copyright 30.11.2011 by Bochkanov Sergey
*************************************************************************/
void linlsqrresults(const linlsqrstate &state, real_1d_array &x, linlsqrrep
ort &rep);
/*************************************************************************
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 MinCGOptimize().
-- ALGLIB --
Copyright 30.11.2011 by Bochkanov Sergey
*************************************************************************/
void linlsqrsetxrep(const linlsqrstate &state, const bool needxrep);
/*************************************************************************
This function initializes linear CG Solver. This solver is used to solve
symmetric positive definite problems. If you want to solve nonsymmetric
(or non-positive definite) problem you may use LinLSQR solver provided by
ALGLIB.
USAGE:
1. User initializes algorithm state with LinCGCreate() call
2. User tunes solver parameters with LinCGSetCond() and other functions
3. Optionally, user sets starting point with LinCGSetStartingPoint()
4. User calls LinCGSolveSparse() function which takes algorithm state and
SparseMatrix object.
5. User calls LinCGResults() to get solution
6. Optionally, user may call LinCGSolveSparse() again to solve another
problem with different matrix and/or right part without reinitializing
LinCGState structure.
INPUT PARAMETERS:
N - problem dimension, N>0
OUTPUT PARAMETERS:
State - structure which stores algorithm state
-- ALGLIB --
Copyright 14.11.2011 by Bochkanov Sergey
*************************************************************************/
void lincgcreate(const ae_int_t n, lincgstate &state);
/*************************************************************************
This function sets starting point.
By default, zero starting point is used.
INPUT PARAMETERS:
X - starting point, array[N]
OUTPUT PARAMETERS:
State - structure which stores algorithm state
-- ALGLIB --
Copyright 14.11.2011 by Bochkanov Sergey
*************************************************************************/
void lincgsetstartingpoint(const lincgstate &state, const real_1d_array &x)
;
/*************************************************************************
This function sets stopping criteria.
INPUT PARAMETERS:
EpsF - algorithm will be stopped if norm of residual is less than
EpsF*||b||.
MaxIts - algorithm will be stopped if number of iterations is more
than MaxIts.
OUTPUT PARAMETERS:
State - structure which stores algorithm state
NOTES:
If both EpsF and MaxIts are zero then small EpsF will be set to small
value.
-- ALGLIB --
Copyright 14.11.2011 by Bochkanov Sergey
*************************************************************************/
void lincgsetcond(const lincgstate &state, const double epsf, const ae_int_
t maxits);
/*************************************************************************
Procedure for solution of A*x=b with sparse A.
INPUT PARAMETERS:
State - algorithm state
A - sparse matrix in the CRS format (you MUST contvert it to
CRS format by calling SparseConvertToCRS() function).
IsUpper - whether upper or lower triangle of A is used:
* IsUpper=True => only upper triangle is used and lower
triangle is not referenced at all
* IsUpper=False => only lower triangle is used and upper
triangle is not referenced at all
B - right part, array[N]
RESULT:
This function returns no result.
You can get solution by calling LinCGResults()
-- ALGLIB --
Copyright 14.11.2011 by Bochkanov Sergey
*************************************************************************/
void lincgsolvesparse(const lincgstate &state, const sparsematrix &a, const
bool isupper, const real_1d_array &b);
/*************************************************************************
CG-solver: results.
This function must be called after LinCGSolve
INPUT PARAMETERS:
State - algorithm state
OUTPUT PARAMETERS:
X - array[N], solution
Rep - optimization report:
* Rep.TerminationType completetion code:
* -5 input matrix is either not positive definite,
too large or too small
* -4 overflow/underflow during solution
(ill conditioned problem)
* 1 ||residual||<=EpsF*||b||
* 5 MaxIts steps was taken
* 7 rounding errors prevent further progress,
best point found is returned
* Rep.IterationsCount contains iterations count
* NMV countains number of matrix-vector calculations
-- ALGLIB --
Copyright 14.11.2011 by Bochkanov Sergey
*************************************************************************/
void lincgresults(const lincgstate &state, real_1d_array &x, lincgreport &r
ep);
/*************************************************************************
This function sets restart frequency. By default, algorithm is restarted
after N subsequent iterations.
-- ALGLIB --
Copyright 14.11.2011 by Bochkanov Sergey
*************************************************************************/
void lincgsetrestartfreq(const lincgstate &state, const ae_int_t srf);
/*************************************************************************
This function sets frequency of residual recalculations.
Algorithm updates residual r_k using iterative formula, but recalculates
it from scratch after each 10 iterations. It is done to avoid accumulation
of numerical errors and to stop algorithm when r_k starts to grow.
Such low update frequence (1/10) gives very little overhead, but makes
algorithm a bit more robust against numerical errors. However, you may
change it
INPUT PARAMETERS:
Freq - desired update frequency, Freq>=0.
Zero value means that no updates will be done.
-- ALGLIB --
Copyright 14.11.2011 by Bochkanov Sergey
*************************************************************************/
void lincgsetrupdatefreq(const lincgstate &state, const ae_int_t freq);
/*************************************************************************
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 MinCGOptimize().
-- ALGLIB --
Copyright 14.11.2011 by Bochkanov Sergey
*************************************************************************/
void lincgsetxrep(const lincgstate &state, const bool needxrep);
/*************************************************************************
LEVENBERG-MARQUARDT-LIKE NONLINEAR SOLVER LEVENBERG-MARQUARDT-LIKE NONLINEAR SOLVER
DESCRIPTION: DESCRIPTION:
This algorithm solves system of nonlinear equations This algorithm solves system of nonlinear equations
F[0](x[0], ..., x[n-1]) = 0 F[0](x[0], ..., x[n-1]) = 0
F[1](x[0], ..., x[n-1]) = 0 F[1](x[0], ..., x[n-1]) = 0
... ...
F[M-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 with M/N do not necessarily coincide. Algorithm converges quadratically
under following conditions: under following conditions:
skipping to change at line 1284 skipping to change at line 1779
ae_int_t* info, ae_int_t* info,
densesolverlsreport* rep, densesolverlsreport* rep,
/* Real */ ae_vector* x, /* Real */ ae_vector* x,
ae_state *_state); ae_state *_state);
ae_bool _densesolverreport_init(densesolverreport* p, ae_state *_state, ae_ bool make_automatic); ae_bool _densesolverreport_init(densesolverreport* p, ae_state *_state, ae_ bool make_automatic);
ae_bool _densesolverreport_init_copy(densesolverreport* dst, densesolverrep ort* src, ae_state *_state, ae_bool make_automatic); ae_bool _densesolverreport_init_copy(densesolverreport* dst, densesolverrep ort* src, ae_state *_state, ae_bool make_automatic);
void _densesolverreport_clear(densesolverreport* p); void _densesolverreport_clear(densesolverreport* p);
ae_bool _densesolverlsreport_init(densesolverlsreport* p, ae_state *_state, ae_bool make_automatic); ae_bool _densesolverlsreport_init(densesolverlsreport* p, ae_state *_state, ae_bool make_automatic);
ae_bool _densesolverlsreport_init_copy(densesolverlsreport* dst, densesolve rlsreport* src, ae_state *_state, ae_bool make_automatic); ae_bool _densesolverlsreport_init_copy(densesolverlsreport* dst, densesolve rlsreport* src, ae_state *_state, ae_bool make_automatic);
void _densesolverlsreport_clear(densesolverlsreport* p); void _densesolverlsreport_clear(densesolverlsreport* p);
void linlsqrcreate(ae_int_t m,
ae_int_t n,
linlsqrstate* state,
ae_state *_state);
void linlsqrsetb(linlsqrstate* state,
/* Real */ ae_vector* b,
ae_state *_state);
void linlsqrsetlambdai(linlsqrstate* state,
double lambdai,
ae_state *_state);
ae_bool linlsqriteration(linlsqrstate* state, ae_state *_state);
void linlsqrsolvesparse(linlsqrstate* state,
sparsematrix* a,
/* Real */ ae_vector* b,
ae_state *_state);
void linlsqrsetcond(linlsqrstate* state,
double epsa,
double epsb,
ae_int_t maxits,
ae_state *_state);
void linlsqrresults(linlsqrstate* state,
/* Real */ ae_vector* x,
linlsqrreport* rep,
ae_state *_state);
void linlsqrsetxrep(linlsqrstate* state,
ae_bool needxrep,
ae_state *_state);
void linlsqrrestart(linlsqrstate* state, ae_state *_state);
ae_bool _linlsqrstate_init(linlsqrstate* p, ae_state *_state, ae_bool make_
automatic);
ae_bool _linlsqrstate_init_copy(linlsqrstate* dst, linlsqrstate* src, ae_st
ate *_state, ae_bool make_automatic);
void _linlsqrstate_clear(linlsqrstate* p);
ae_bool _linlsqrreport_init(linlsqrreport* p, ae_state *_state, ae_bool mak
e_automatic);
ae_bool _linlsqrreport_init_copy(linlsqrreport* dst, linlsqrreport* src, ae
_state *_state, ae_bool make_automatic);
void _linlsqrreport_clear(linlsqrreport* p);
void lincgcreate(ae_int_t n, lincgstate* state, ae_state *_state);
void lincgsetstartingpoint(lincgstate* state,
/* Real */ ae_vector* x,
ae_state *_state);
void lincgsetb(lincgstate* state,
/* Real */ ae_vector* b,
ae_state *_state);
void lincgsetcond(lincgstate* state,
double epsf,
ae_int_t maxits,
ae_state *_state);
ae_bool lincgiteration(lincgstate* state, ae_state *_state);
void lincgsolvesparse(lincgstate* state,
sparsematrix* a,
ae_bool isupper,
/* Real */ ae_vector* b,
ae_state *_state);
void lincgresults(lincgstate* state,
/* Real */ ae_vector* x,
lincgreport* rep,
ae_state *_state);
void lincgsetrestartfreq(lincgstate* state,
ae_int_t srf,
ae_state *_state);
void lincgsetrupdatefreq(lincgstate* state,
ae_int_t freq,
ae_state *_state);
void lincgsetxrep(lincgstate* state, ae_bool needxrep, ae_state *_state);
void lincgrestart(lincgstate* state, ae_state *_state);
ae_bool _lincgstate_init(lincgstate* p, ae_state *_state, ae_bool make_auto
matic);
ae_bool _lincgstate_init_copy(lincgstate* dst, lincgstate* src, ae_state *_
state, ae_bool make_automatic);
void _lincgstate_clear(lincgstate* p);
ae_bool _lincgreport_init(lincgreport* p, ae_state *_state, ae_bool make_au
tomatic);
ae_bool _lincgreport_init_copy(lincgreport* dst, lincgreport* src, ae_state
*_state, ae_bool make_automatic);
void _lincgreport_clear(lincgreport* p);
void nleqcreatelm(ae_int_t n, void nleqcreatelm(ae_int_t n,
ae_int_t m, ae_int_t m,
/* Real */ ae_vector* x, /* Real */ ae_vector* x,
nleqstate* state, nleqstate* state,
ae_state *_state); ae_state *_state);
void nleqsetcond(nleqstate* state, void nleqsetcond(nleqstate* state,
double epsf, double epsf,
ae_int_t maxits, ae_int_t maxits,
ae_state *_state); ae_state *_state);
void nleqsetxrep(nleqstate* state, ae_bool needxrep, ae_state *_state); void nleqsetxrep(nleqstate* state, ae_bool needxrep, ae_state *_state);
 End of changes. 4 change blocks. 
0 lines changed or deleted 581 lines changed or added


 statistics.h   statistics.h 
skipping to change at line 66 skipping to change at line 66
Skewness- skewness (if variance<>0; zero otherwise). Skewness- skewness (if variance<>0; zero otherwise).
Kurtosis- kurtosis (if variance<>0; zero otherwise). Kurtosis- kurtosis (if variance<>0; zero otherwise).
-- ALGLIB -- -- ALGLIB --
Copyright 06.09.2006 by Bochkanov Sergey Copyright 06.09.2006 by Bochkanov Sergey
*************************************************************************/ *************************************************************************/
void samplemoments(const real_1d_array &x, const ae_int_t n, double &mean, double &variance, double &skewness, double &kurtosis); void samplemoments(const real_1d_array &x, const ae_int_t n, double &mean, double &variance, double &skewness, double &kurtosis);
void samplemoments(const real_1d_array &x, double &mean, double &variance, double &skewness, double &kurtosis); void samplemoments(const real_1d_array &x, double &mean, double &variance, double &skewness, double &kurtosis);
/************************************************************************* /*************************************************************************
Calculation of the mean.
INPUT PARAMETERS:
X - sample
N - N>=0, sample size:
* if given, only leading N elements of X are processed
* if not given, automatically determined from size of X
NOTE:
This function return result which calculated by 'SampleMoments' function
and stored at 'Mean' variable.
-- ALGLIB --
Copyright 06.09.2006 by Bochkanov Sergey
*************************************************************************/
double samplemean(const real_1d_array &x, const ae_int_t n);
double samplemean(const real_1d_array &x);
/*************************************************************************
Calculation of the variance.
INPUT PARAMETERS:
X - sample
N - N>=0, sample size:
* if given, only leading N elements of X are processed
* if not given, automatically determined from size of X
NOTE:
This function return result which calculated by 'SampleMoments' function
and stored at 'Variance' variable.
-- ALGLIB --
Copyright 06.09.2006 by Bochkanov Sergey
*************************************************************************/
double samplevariance(const real_1d_array &x, const ae_int_t n);
double samplevariance(const real_1d_array &x);
/*************************************************************************
Calculation of the skewness.
INPUT PARAMETERS:
X - sample
N - N>=0, sample size:
* if given, only leading N elements of X are processed
* if not given, automatically determined from size of X
NOTE:
This function return result which calculated by 'SampleMoments' function
and stored at 'Skewness' variable.
-- ALGLIB --
Copyright 06.09.2006 by Bochkanov Sergey
*************************************************************************/
double sampleskewness(const real_1d_array &x, const ae_int_t n);
double sampleskewness(const real_1d_array &x);
/*************************************************************************
Calculation of the kurtosis.
INPUT PARAMETERS:
X - sample
N - N>=0, sample size:
* if given, only leading N elements of X are processed
* if not given, automatically determined from size of X
NOTE:
This function return result which calculated by 'SampleMoments' function
and stored at 'Kurtosis' variable.
-- ALGLIB --
Copyright 06.09.2006 by Bochkanov Sergey
*************************************************************************/
double samplekurtosis(const real_1d_array &x, const ae_int_t n);
double samplekurtosis(const real_1d_array &x);
/*************************************************************************
ADev ADev
Input parameters: Input parameters:
X - sample X - sample
N - N>=0, sample size: N - N>=0, sample size:
* if given, only leading N elements of X are processed * if given, only leading N elements of X are processed
* if not given, automatically determined from size of X * if not given, automatically determined from size of X
Output parameters: Output parameters:
ADev- ADev ADev- ADev
skipping to change at line 830 skipping to change at line 910
///////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////
namespace alglib_impl namespace alglib_impl
{ {
void samplemoments(/* Real */ ae_vector* x, void samplemoments(/* Real */ ae_vector* x,
ae_int_t n, ae_int_t n,
double* mean, double* mean,
double* variance, double* variance,
double* skewness, double* skewness,
double* kurtosis, double* kurtosis,
ae_state *_state); ae_state *_state);
double samplemean(/* Real */ ae_vector* x,
ae_int_t n,
ae_state *_state);
double samplevariance(/* Real */ ae_vector* x,
ae_int_t n,
ae_state *_state);
double sampleskewness(/* Real */ ae_vector* x,
ae_int_t n,
ae_state *_state);
double samplekurtosis(/* Real */ ae_vector* x,
ae_int_t n,
ae_state *_state);
void sampleadev(/* Real */ ae_vector* x, void sampleadev(/* Real */ ae_vector* x,
ae_int_t n, ae_int_t n,
double* adev, double* adev,
ae_state *_state); ae_state *_state);
void samplemedian(/* Real */ ae_vector* x, void samplemedian(/* Real */ ae_vector* x,
ae_int_t n, ae_int_t n,
double* median, double* median,
ae_state *_state); ae_state *_state);
void samplepercentile(/* Real */ ae_vector* x, void samplepercentile(/* Real */ ae_vector* x,
ae_int_t n, ae_int_t n,
 End of changes. 2 change blocks. 
0 lines changed or deleted 92 lines changed or added

This html diff was produced by rfcdiff 1.41. The latest version is available from http://tools.ietf.org/tools/rfcdiff/