| interpolation.h | | interpolation.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 _interpolation_pkg_h | | #ifndef _interpolation_pkg_h | |
| #define _interpolation_pkg_h | | #define _interpolation_pkg_h | |
| #include "ap.h" | | #include "ap.h" | |
| #include "alglibinternal.h" | | #include "alglibinternal.h" | |
| #include "alglibmisc.h" | | #include "alglibmisc.h" | |
| #include "linalg.h" | | #include "linalg.h" | |
|
| #include "solvers.h" | | | |
| #include "optimization.h" | | #include "optimization.h" | |
|
| | | #include "solvers.h" | |
| #include "specialfunctions.h" | | #include "specialfunctions.h" | |
| #include "integration.h" | | #include "integration.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 | |
| { | | { | |
| double taskrcond; | | double taskrcond; | |
| double rmserror; | | double rmserror; | |
| double avgerror; | | double avgerror; | |
| double avgrelerror; | | double avgrelerror; | |
| double maxerror; | | double maxerror; | |
|
| | | } polynomialfitreport; | |
| | | typedef struct | |
| | | { | |
| | | double taskrcond; | |
| | | ae_int_t dbest; | |
| | | double rmserror; | |
| | | double avgerror; | |
| | | double avgrelerror; | |
| | | double maxerror; | |
| | | } barycentricfitreport; | |
| | | typedef struct | |
| | | { | |
| | | double taskrcond; | |
| | | double rmserror; | |
| | | double avgerror; | |
| | | double avgrelerror; | |
| | | double maxerror; | |
| | | } spline1dfitreport; | |
| | | typedef struct | |
| | | { | |
| | | double taskrcond; | |
| | | double rmserror; | |
| | | double avgerror; | |
| | | double avgrelerror; | |
| | | double maxerror; | |
| } lsfitreport; | | } lsfitreport; | |
| typedef struct | | typedef struct | |
| { | | { | |
| ae_int_t n; | | ae_int_t n; | |
| ae_int_t m; | | ae_int_t m; | |
| ae_int_t k; | | ae_int_t k; | |
| double epsf; | | double epsf; | |
| double epsx; | | double epsx; | |
| ae_int_t maxits; | | ae_int_t maxits; | |
| double stpmax; | | double stpmax; | |
| ae_bool xrep; | | ae_bool xrep; | |
| ae_matrix taskx; | | ae_matrix taskx; | |
| ae_vector tasky; | | ae_vector tasky; | |
| ae_vector w; | | ae_vector w; | |
|
| ae_bool cheapfg; | | | |
| ae_bool havehess; | | | |
| ae_bool xupdated; | | ae_bool xupdated; | |
| ae_bool needf; | | ae_bool needf; | |
| ae_bool needfg; | | ae_bool needfg; | |
| ae_bool needfgh; | | ae_bool needfgh; | |
| ae_int_t pointindex; | | ae_int_t pointindex; | |
| ae_vector x; | | ae_vector x; | |
| ae_vector c; | | ae_vector c; | |
| double f; | | double f; | |
| ae_vector g; | | ae_vector g; | |
| ae_matrix h; | | ae_matrix h; | |
| | | | |
| skipping to change at line 81 | | skipping to change at line 104 | |
|---|
| double reprmserror; | | double reprmserror; | |
| double repavgerror; | | double repavgerror; | |
| double repavgrelerror; | | double repavgrelerror; | |
| double repmaxerror; | | double repmaxerror; | |
| minlmstate optstate; | | minlmstate optstate; | |
| minlmreport optrep; | | minlmreport optrep; | |
| rcommstate rstate; | | rcommstate rstate; | |
| } lsfitstate; | | } lsfitstate; | |
| typedef struct | | typedef struct | |
| { | | { | |
|
| double taskrcond; | | | |
| double rmserror; | | | |
| double avgerror; | | | |
| double avgrelerror; | | | |
| double maxerror; | | | |
| } polynomialfitreport; | | | |
| typedef struct | | | |
| { | | | |
| ae_int_t n; | | ae_int_t n; | |
| 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 | |
| { | | { | |
|
| double taskrcond; | | | |
| ae_int_t dbest; | | | |
| double rmserror; | | | |
| double avgerror; | | | |
| double avgrelerror; | | | |
| double maxerror; | | | |
| } barycentricfitreport; | | | |
| typedef struct | | | |
| { | | | |
| ae_int_t k; | | | |
| ae_vector c; | | | |
| } spline2dinterpolant; | | | |
| 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_vector x; | | ae_vector x; | |
| ae_vector c; | | ae_vector c; | |
| } spline1dinterpolant; | | } spline1dinterpolant; | |
| typedef struct | | typedef struct | |
| { | | { | |
|
| double taskrcond; | | ae_int_t k; | |
| double rmserror; | | ae_vector c; | |
| double avgerror; | | } spline2dinterpolant; | |
| double avgrelerror; | | | |
| double maxerror; | | | |
| } spline1dfitreport; | | | |
| typedef struct | | typedef struct | |
| { | | { | |
| ae_int_t n; | | ae_int_t n; | |
| ae_int_t nx; | | ae_int_t nx; | |
| ae_int_t d; | | ae_int_t d; | |
| double r; | | double r; | |
| ae_int_t nw; | | ae_int_t nw; | |
| kdtree tree; | | kdtree tree; | |
| ae_int_t modeltype; | | ae_int_t modeltype; | |
| ae_matrix q; | | ae_matrix q; | |
| | | | |
| skipping to change at line 172 | | skipping to change at line 170 | |
|---|
| | | | |
| ///////////////////////////////////////////////////////////////////////// | | ///////////////////////////////////////////////////////////////////////// | |
| // | | // | |
| // THIS SECTION CONTAINS C++ INTERFACE | | // THIS SECTION CONTAINS C++ INTERFACE | |
| // | | // | |
| ///////////////////////////////////////////////////////////////////////// | | ///////////////////////////////////////////////////////////////////////// | |
| namespace alglib | | namespace alglib | |
| { | | { | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
|
| | | Polynomial fitting report: | |
| | | TaskRCond reciprocal of task's condition number | |
| | | RMSError RMS error | |
| | | AvgError average error | |
| | | AvgRelError average relative error (for non-zero Y[I]) | |
| | | MaxError maximum error | |
| | | *************************************************************************/ | |
| | | class _polynomialfitreport_owner | |
| | | { | |
| | | public: | |
| | | _polynomialfitreport_owner(); | |
| | | _polynomialfitreport_owner(const _polynomialfitreport_owner &rhs); | |
| | | _polynomialfitreport_owner& operator=(const _polynomialfitreport_owner | |
| | | &rhs); | |
| | | virtual ~_polynomialfitreport_owner(); | |
| | | alglib_impl::polynomialfitreport* c_ptr(); | |
| | | alglib_impl::polynomialfitreport* c_ptr() const; | |
| | | protected: | |
| | | alglib_impl::polynomialfitreport *p_struct; | |
| | | }; | |
| | | class polynomialfitreport : public _polynomialfitreport_owner | |
| | | { | |
| | | public: | |
| | | polynomialfitreport(); | |
| | | polynomialfitreport(const polynomialfitreport &rhs); | |
| | | polynomialfitreport& operator=(const polynomialfitreport &rhs); | |
| | | virtual ~polynomialfitreport(); | |
| | | double &taskrcond; | |
| | | double &rmserror; | |
| | | double &avgerror; | |
| | | double &avgrelerror; | |
| | | double &maxerror; | |
| | | | |
| | | }; | |
| | | | |
| | | /************************************************************************* | |
| | | Barycentric fitting report: | |
| | | RMSError RMS error | |
| | | AvgError average error | |
| | | AvgRelError average relative error (for non-zero Y[I]) | |
| | | MaxError maximum error | |
| | | TaskRCond reciprocal of task's condition number | |
| | | *************************************************************************/ | |
| | | class _barycentricfitreport_owner | |
| | | { | |
| | | public: | |
| | | _barycentricfitreport_owner(); | |
| | | _barycentricfitreport_owner(const _barycentricfitreport_owner &rhs); | |
| | | _barycentricfitreport_owner& operator=(const _barycentricfitreport_owne | |
| | | r &rhs); | |
| | | virtual ~_barycentricfitreport_owner(); | |
| | | alglib_impl::barycentricfitreport* c_ptr(); | |
| | | alglib_impl::barycentricfitreport* c_ptr() const; | |
| | | protected: | |
| | | alglib_impl::barycentricfitreport *p_struct; | |
| | | }; | |
| | | class barycentricfitreport : public _barycentricfitreport_owner | |
| | | { | |
| | | public: | |
| | | barycentricfitreport(); | |
| | | barycentricfitreport(const barycentricfitreport &rhs); | |
| | | barycentricfitreport& operator=(const barycentricfitreport &rhs); | |
| | | virtual ~barycentricfitreport(); | |
| | | double &taskrcond; | |
| | | ae_int_t &dbest; | |
| | | double &rmserror; | |
| | | double &avgerror; | |
| | | double &avgrelerror; | |
| | | double &maxerror; | |
| | | | |
| | | }; | |
| | | | |
| | | /************************************************************************* | |
| | | Spline fitting report: | |
| | | RMSError RMS error | |
| | | AvgError average error | |
| | | AvgRelError average relative error (for non-zero Y[I]) | |
| | | MaxError maximum error | |
| | | | |
| | | Fields below are filled by obsolete functions (Spline1DFitCubic, | |
| | | Spline1DFitHermite). Modern fitting functions do NOT fill these fields: | |
| | | TaskRCond reciprocal of task's condition number | |
| | | *************************************************************************/ | |
| | | class _spline1dfitreport_owner | |
| | | { | |
| | | public: | |
| | | _spline1dfitreport_owner(); | |
| | | _spline1dfitreport_owner(const _spline1dfitreport_owner &rhs); | |
| | | _spline1dfitreport_owner& operator=(const _spline1dfitreport_owner &rhs | |
| | | ); | |
| | | virtual ~_spline1dfitreport_owner(); | |
| | | alglib_impl::spline1dfitreport* c_ptr(); | |
| | | alglib_impl::spline1dfitreport* c_ptr() const; | |
| | | protected: | |
| | | alglib_impl::spline1dfitreport *p_struct; | |
| | | }; | |
| | | class spline1dfitreport : public _spline1dfitreport_owner | |
| | | { | |
| | | public: | |
| | | spline1dfitreport(); | |
| | | spline1dfitreport(const spline1dfitreport &rhs); | |
| | | spline1dfitreport& operator=(const spline1dfitreport &rhs); | |
| | | virtual ~spline1dfitreport(); | |
| | | double &taskrcond; | |
| | | double &rmserror; | |
| | | double &avgerror; | |
| | | double &avgrelerror; | |
| | | double &maxerror; | |
| | | | |
| | | }; | |
| | | | |
| | | /************************************************************************* | |
| Least squares fitting report: | | Least squares fitting report: | |
| TaskRCond reciprocal of task's condition number | | TaskRCond reciprocal of task's condition number | |
| RMSError RMS error | | RMSError RMS error | |
| AvgError average error | | AvgError average error | |
| AvgRelError average relative error (for non-zero Y[I]) | | AvgRelError average relative error (for non-zero Y[I]) | |
| MaxError maximum error | | MaxError maximum error | |
| *************************************************************************/ | | *************************************************************************/ | |
| class _lsfitreport_owner | | class _lsfitreport_owner | |
| { | | { | |
| public: | | public: | |
| | | | |
| skipping to change at line 207 | | skipping to change at line 314 | |
|---|
| virtual ~lsfitreport(); | | virtual ~lsfitreport(); | |
| double &taskrcond; | | double &taskrcond; | |
| double &rmserror; | | double &rmserror; | |
| double &avgerror; | | double &avgerror; | |
| double &avgrelerror; | | double &avgrelerror; | |
| double &maxerror; | | double &maxerror; | |
| | | | |
| }; | | }; | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
|
| | | Nonlinear fitter. | |
| | | | |
|
| | | You should use ALGLIB functions to work with fitter. | |
| | | Never try to access its fields directly! | |
| *************************************************************************/ | | *************************************************************************/ | |
| class _lsfitstate_owner | | class _lsfitstate_owner | |
| { | | { | |
| public: | | public: | |
| _lsfitstate_owner(); | | _lsfitstate_owner(); | |
| _lsfitstate_owner(const _lsfitstate_owner &rhs); | | _lsfitstate_owner(const _lsfitstate_owner &rhs); | |
| _lsfitstate_owner& operator=(const _lsfitstate_owner &rhs); | | _lsfitstate_owner& operator=(const _lsfitstate_owner &rhs); | |
| virtual ~_lsfitstate_owner(); | | virtual ~_lsfitstate_owner(); | |
| alglib_impl::lsfitstate* c_ptr(); | | alglib_impl::lsfitstate* c_ptr(); | |
| alglib_impl::lsfitstate* c_ptr() const; | | alglib_impl::lsfitstate* c_ptr() const; | |
| | | | |
| skipping to change at line 241 | | skipping to change at line 351 | |
|---|
| ae_bool &xupdated; | | ae_bool &xupdated; | |
| real_1d_array c; | | real_1d_array c; | |
| double &f; | | double &f; | |
| real_1d_array g; | | real_1d_array g; | |
| real_2d_array h; | | real_2d_array h; | |
| real_1d_array x; | | real_1d_array x; | |
| | | | |
| }; | | }; | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
|
| Polynomial fitting report: | | | |
| TaskRCond reciprocal of task's condition number | | | |
| RMSError RMS error | | | |
| AvgError average error | | | |
| AvgRelError average relative error (for non-zero Y[I]) | | | |
| MaxError maximum error | | | |
| *************************************************************************/ | | | |
| class _polynomialfitreport_owner | | | |
| { | | | |
| public: | | | |
| _polynomialfitreport_owner(); | | | |
| _polynomialfitreport_owner(const _polynomialfitreport_owner &rhs); | | | |
| _polynomialfitreport_owner& operator=(const _polynomialfitreport_owner | | | |
| &rhs); | | | |
| virtual ~_polynomialfitreport_owner(); | | | |
| alglib_impl::polynomialfitreport* c_ptr(); | | | |
| alglib_impl::polynomialfitreport* c_ptr() const; | | | |
| protected: | | | |
| alglib_impl::polynomialfitreport *p_struct; | | | |
| }; | | | |
| class polynomialfitreport : public _polynomialfitreport_owner | | | |
| { | | | |
| public: | | | |
| polynomialfitreport(); | | | |
| polynomialfitreport(const polynomialfitreport &rhs); | | | |
| polynomialfitreport& operator=(const polynomialfitreport &rhs); | | | |
| virtual ~polynomialfitreport(); | | | |
| double &taskrcond; | | | |
| double &rmserror; | | | |
| double &avgerror; | | | |
| double &avgrelerror; | | | |
| double &maxerror; | | | |
| | | | |
| }; | | | |
| | | | |
| /************************************************************************* | | | |
| Barycentric interpolant. | | Barycentric interpolant. | |
| *************************************************************************/ | | *************************************************************************/ | |
| class _barycentricinterpolant_owner | | class _barycentricinterpolant_owner | |
| { | | { | |
| public: | | public: | |
| _barycentricinterpolant_owner(); | | _barycentricinterpolant_owner(); | |
| _barycentricinterpolant_owner(const _barycentricinterpolant_owner &rhs)
; | | _barycentricinterpolant_owner(const _barycentricinterpolant_owner &rhs)
; | |
| _barycentricinterpolant_owner& operator=(const _barycentricinterpolant_
owner &rhs); | | _barycentricinterpolant_owner& operator=(const _barycentricinterpolant_
owner &rhs); | |
| virtual ~_barycentricinterpolant_owner(); | | virtual ~_barycentricinterpolant_owner(); | |
| alglib_impl::barycentricinterpolant* c_ptr(); | | alglib_impl::barycentricinterpolant* c_ptr(); | |
| | | | |
| skipping to change at line 301 | | skipping to change at line 376 | |
|---|
| { | | { | |
| 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(); | |
| | | | |
| }; | | }; | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
|
| Barycentric fitting report: | | 1-dimensional spline inteprolant | |
| TaskRCond reciprocal of task's condition number | | | |
| RMSError RMS error | | | |
| AvgError average error | | | |
| AvgRelError average relative error (for non-zero Y[I]) | | | |
| MaxError maximum error | | | |
| *************************************************************************/ | | *************************************************************************/ | |
|
| class _barycentricfitreport_owner | | class _spline1dinterpolant_owner | |
| { | | { | |
| public: | | public: | |
|
| _barycentricfitreport_owner(); | | _spline1dinterpolant_owner(); | |
| _barycentricfitreport_owner(const _barycentricfitreport_owner &rhs); | | _spline1dinterpolant_owner(const _spline1dinterpolant_owner &rhs); | |
| _barycentricfitreport_owner& operator=(const _barycentricfitreport_owne | | _spline1dinterpolant_owner& operator=(const _spline1dinterpolant_owner | |
| r &rhs); | | &rhs); | |
| virtual ~_barycentricfitreport_owner(); | | virtual ~_spline1dinterpolant_owner(); | |
| alglib_impl::barycentricfitreport* c_ptr(); | | alglib_impl::spline1dinterpolant* c_ptr(); | |
| alglib_impl::barycentricfitreport* c_ptr() const; | | alglib_impl::spline1dinterpolant* c_ptr() const; | |
| protected: | | protected: | |
|
| alglib_impl::barycentricfitreport *p_struct; | | alglib_impl::spline1dinterpolant *p_struct; | |
| }; | | }; | |
|
| class barycentricfitreport : public _barycentricfitreport_owner | | class spline1dinterpolant : public _spline1dinterpolant_owner | |
| { | | { | |
| public: | | public: | |
|
| barycentricfitreport(); | | spline1dinterpolant(); | |
| barycentricfitreport(const barycentricfitreport &rhs); | | spline1dinterpolant(const spline1dinterpolant &rhs); | |
| barycentricfitreport& operator=(const barycentricfitreport &rhs); | | spline1dinterpolant& operator=(const spline1dinterpolant &rhs); | |
| virtual ~barycentricfitreport(); | | virtual ~spline1dinterpolant(); | |
| double &taskrcond; | | | |
| ae_int_t &dbest; | | | |
| double &rmserror; | | | |
| double &avgerror; | | | |
| double &avgrelerror; | | | |
| double &maxerror; | | | |
| | | | |
| }; | | }; | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
| 2-dimensional spline inteprolant | | 2-dimensional spline inteprolant | |
| *************************************************************************/ | | *************************************************************************/ | |
| class _spline2dinterpolant_owner | | class _spline2dinterpolant_owner | |
| { | | { | |
| public: | | public: | |
| _spline2dinterpolant_owner(); | | _spline2dinterpolant_owner(); | |
| | | | |
| skipping to change at line 362 | | skipping to change at line 426 | |
|---|
| { | | { | |
| public: | | public: | |
| spline2dinterpolant(); | | spline2dinterpolant(); | |
| spline2dinterpolant(const spline2dinterpolant &rhs); | | spline2dinterpolant(const spline2dinterpolant &rhs); | |
| spline2dinterpolant& operator=(const spline2dinterpolant &rhs); | | spline2dinterpolant& operator=(const spline2dinterpolant &rhs); | |
| virtual ~spline2dinterpolant(); | | virtual ~spline2dinterpolant(); | |
| | | | |
| }; | | }; | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
|
| 1-dimensional spline inteprolant | | IDW interpolant. | |
| *************************************************************************/ | | *************************************************************************/ | |
|
| class _spline1dinterpolant_owner | | class _idwinterpolant_owner | |
| { | | { | |
| public: | | public: | |
|
| _spline1dinterpolant_owner(); | | _idwinterpolant_owner(); | |
| _spline1dinterpolant_owner(const _spline1dinterpolant_owner &rhs); | | _idwinterpolant_owner(const _idwinterpolant_owner &rhs); | |
| _spline1dinterpolant_owner& operator=(const _spline1dinterpolant_owner | | _idwinterpolant_owner& operator=(const _idwinterpolant_owner &rhs); | |
| &rhs); | | virtual ~_idwinterpolant_owner(); | |
| virtual ~_spline1dinterpolant_owner(); | | alglib_impl::idwinterpolant* c_ptr(); | |
| alglib_impl::spline1dinterpolant* c_ptr(); | | alglib_impl::idwinterpolant* c_ptr() const; | |
| alglib_impl::spline1dinterpolant* c_ptr() const; | | | |
| protected: | | | |
| alglib_impl::spline1dinterpolant *p_struct; | | | |
| }; | | | |
| class spline1dinterpolant : public _spline1dinterpolant_owner | | | |
| { | | | |
| public: | | | |
| spline1dinterpolant(); | | | |
| spline1dinterpolant(const spline1dinterpolant &rhs); | | | |
| spline1dinterpolant& operator=(const spline1dinterpolant &rhs); | | | |
| virtual ~spline1dinterpolant(); | | | |
| | | | |
| }; | | | |
| | | | |
| /************************************************************************* | | | |
| Spline fitting report: | | | |
| TaskRCond reciprocal of task's condition number | | | |
| RMSError RMS error | | | |
| AvgError average error | | | |
| AvgRelError average relative error (for non-zero Y[I]) | | | |
| MaxError maximum error | | | |
| *************************************************************************/ | | | |
| class _spline1dfitreport_owner | | | |
| { | | | |
| public: | | | |
| _spline1dfitreport_owner(); | | | |
| _spline1dfitreport_owner(const _spline1dfitreport_owner &rhs); | | | |
| _spline1dfitreport_owner& operator=(const _spline1dfitreport_owner &rhs | | | |
| ); | | | |
| virtual ~_spline1dfitreport_owner(); | | | |
| alglib_impl::spline1dfitreport* c_ptr(); | | | |
| alglib_impl::spline1dfitreport* c_ptr() const; | | | |
| protected: | | | |
| alglib_impl::spline1dfitreport *p_struct; | | | |
| }; | | | |
| class spline1dfitreport : public _spline1dfitreport_owner | | | |
| { | | | |
| public: | | | |
| spline1dfitreport(); | | | |
| spline1dfitreport(const spline1dfitreport &rhs); | | | |
| spline1dfitreport& operator=(const spline1dfitreport &rhs); | | | |
| virtual ~spline1dfitreport(); | | | |
| double &taskrcond; | | | |
| double &rmserror; | | | |
| double &avgerror; | | | |
| double &avgrelerror; | | | |
| double &maxerror; | | | |
| | | | |
| }; | | | |
| | | | |
| /************************************************************************* | | | |
| IDW interpolant. | | | |
| *************************************************************************/ | | | |
| class _idwinterpolant_owner | | | |
| { | | | |
| public: | | | |
| _idwinterpolant_owner(); | | | |
| _idwinterpolant_owner(const _idwinterpolant_owner &rhs); | | | |
| _idwinterpolant_owner& operator=(const _idwinterpolant_owner &rhs); | | | |
| virtual ~_idwinterpolant_owner(); | | | |
| alglib_impl::idwinterpolant* c_ptr(); | | | |
| alglib_impl::idwinterpolant* c_ptr() const; | | | |
| protected: | | protected: | |
| alglib_impl::idwinterpolant *p_struct; | | alglib_impl::idwinterpolant *p_struct; | |
| }; | | }; | |
| class idwinterpolant : public _idwinterpolant_owner | | class idwinterpolant : public _idwinterpolant_owner | |
| { | | { | |
| public: | | public: | |
| idwinterpolant(); | | idwinterpolant(); | |
| idwinterpolant(const idwinterpolant &rhs); | | idwinterpolant(const idwinterpolant &rhs); | |
| idwinterpolant& operator=(const idwinterpolant &rhs); | | idwinterpolant& operator=(const idwinterpolant &rhs); | |
| virtual ~idwinterpolant(); | | virtual ~idwinterpolant(); | |
| | | | |
| skipping to change at line 503 | | skipping to change at line 507 | |
|---|
| { | | { | |
| 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(); | |
| | | | |
| }; | | }; | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
|
| Weighted linear least squares fitting. | | Fitting by polynomials in barycentric form. This function provides simple | |
| | | unterface for unconstrained unweighted fitting. See PolynomialFitWC() if | |
| | | you need constrained fitting. | |
| | | | |
|
| QR decomposition is used to reduce task to MxM, then triangular solver or | | Task is linear, so linear least squares solver is used. Complexity of this | |
| SVD-based solver is used depending on condition number of the system. It | | computational scheme is O(N*M^2), mostly dominated by least squares solver | |
| allows to maximize speed and retain decent accuracy. | | | |
| | | SEE ALSO: | |
| | | PolynomialFitWC() | |
| | | | |
| INPUT PARAMETERS: | | INPUT PARAMETERS: | |
|
| Y - array[0..N-1] Function values in N points. | | X - points, array[0..N-1]. | |
| W - array[0..N-1] Weights corresponding to function values. | | Y - function values, array[0..N-1]. | |
| Each summand in square sum of approximation deviations | | N - number of points, N>0 | |
| from given values is multiplied by the square of | | * if given, only leading N elements of X/Y are used | |
| corresponding weight. | | * if not given, automatically determined from sizes of X/Y | |
| FMatrix - a table of basis functions values, array[0..N-1, 0..M-1]. | | M - number of basis functions (= polynomial_degree + 1), M>=1 | |
| FMatrix[I, J] - value of J-th basis function in I-th point. | | | |
| N - number of points used. N>=1. | | | |
| M - number of basis functions, M>=1. | | | |
| | | | |
| OUTPUT PARAMETERS: | | OUTPUT PARAMETERS: | |
|
| Info - error code: | | Info- same format as in LSFitLinearW() subroutine: | |
| * -4 internal SVD decomposition subroutine failed (very | | * Info>0 task is solved | |
| rare and for degenerate systems only) | | * Info<=0 an error occured: | |
| * -1 incorrect N/M were specified | | -4 means inconvergence of internal SVD | |
| * 1 task is solved | | P - interpolant in barycentric form. | |
| C - decomposition coefficients, array[0..M-1] | | Rep - report, same format as in LSFitLinearW() subroutine. | |
| Rep - fitting report. Following fields are set: | | Following fields are set: | |
| * Rep.TaskRCond reciprocal of condition number | | * RMSError rms error on the (X,Y). | |
| * RMSError rms error on the (X,Y). | | * AvgError average error on the (X,Y). | |
| * AvgError average error on the (X,Y). | | * AvgRelError average relative error on the non-zero Y | |
| * AvgRelError average relative error on the non-zero | | * MaxError maximum error | |
| Y | | NON-WEIGHTED ERRORS ARE CALCULATED | |
| * MaxError maximum error | | | |
| NON-WEIGHTED ERRORS ARE CALCULATED | | | |
| | | | |
|
| SEE ALSO | | NOTES: | |
| LSFitLinear | | you can convert P from barycentric form to the power or Chebyshev | |
| LSFitLinearC | | basis with PolynomialBar2Pow() or PolynomialBar2Cheb() functions from | |
| LSFitLinearWC | | POLINT subpackage. | |
| | | | |
|
| -- ALGLIB -- | | -- ALGLIB PROJECT -- | |
| Copyright 17.08.2009 by Bochkanov Sergey | | Copyright 10.12.2009 by Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
|
| void lsfitlinearw(const real_1d_array &y, const real_1d_array &w, const rea | | void polynomialfit(const real_1d_array &x, const real_1d_array &y, const ae | |
| l_2d_array &fmatrix, const ae_int_t n, const ae_int_t m, ae_int_t &info, re | | _int_t n, const ae_int_t m, ae_int_t &info, barycentricinterpolant &p, poly | |
| al_1d_array &c, lsfitreport &rep); | | nomialfitreport &rep); | |
| void lsfitlinearw(const real_1d_array &y, const real_1d_array &w, const rea | | void polynomialfit(const real_1d_array &x, const real_1d_array &y, const ae | |
| l_2d_array &fmatrix, ae_int_t &info, real_1d_array &c, lsfitreport &rep); | | _int_t m, ae_int_t &info, barycentricinterpolant &p, polynomialfitreport &r | |
| | | ep); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
|
| Weighted constained linear least squares fitting. | | Weighted fitting by polynomials in barycentric form, with constraints on | |
| | | function values or first derivatives. | |
| | | | |
|
| This is variation of LSFitLinearW(), which searchs for min|A*x=b| given | | Small regularizing term is used when solving constrained tasks (to improve | |
| that K additional constaints C*x=bc are satisfied. It reduces original | | stability). | |
| task to modified one: min|B*y-d| WITHOUT constraints, then LSFitLinearW() | | | |
| is called. | | Task is linear, so linear least squares solver is used. Complexity of this | |
| | | computational scheme is O(N*M^2), mostly dominated by least squares solver | |
| | | | |
| | | SEE ALSO: | |
| | | PolynomialFit() | |
| | | | |
| INPUT PARAMETERS: | | INPUT PARAMETERS: | |
|
| Y - array[0..N-1] Function values in N points. | | X - points, array[0..N-1]. | |
| W - array[0..N-1] Weights corresponding to function values. | | Y - function values, array[0..N-1]. | |
| Each summand in square sum of approximation deviations | | W - weights, array[0..N-1] | |
| from given values is multiplied by the square of | | Each summand in square sum of approximation deviations from | |
| corresponding weight. | | given values is multiplied by the square of corresponding | |
| FMatrix - a table of basis functions values, array[0..N-1, 0..M-1]. | | weight. Fill it by 1's if you don't want to solve weighted | |
| FMatrix[I,J] - value of J-th basis function in I-th point. | | task. | |
| CMatrix - a table of constaints, array[0..K-1,0..M]. | | N - number of points, N>0. | |
| I-th row of CMatrix corresponds to I-th linear constraint: | | * if given, only leading N elements of X/Y/W are used | |
| CMatrix[I,0]*C[0] + ... + CMatrix[I,M-1]*C[M-1] = CMatrix[I | | * if not given, automatically determined from sizes of X/Y/W | |
| ,M] | | XC - points where polynomial values/derivatives are constrained, | |
| N - number of points used. N>=1. | | array[0..K-1]. | |
| M - number of basis functions, M>=1. | | YC - values of constraints, array[0..K-1] | |
| K - number of constraints, 0 <= K < M | | DC - array[0..K-1], types of constraints: | |
| K=0 corresponds to absence of constraints. | | * DC[i]=0 means that P(XC[i])=YC[i] | |
| | | * DC[i]=1 means that P'(XC[i])=YC[i] | |
| | | SEE BELOW FOR IMPORTANT INFORMATION ON CONSTRAINTS | |
| | | K - number of constraints, 0<=K<M. | |
| | | K=0 means no constraints (XC/YC/DC are not used in such cases) | |
| | | M - number of basis functions (= polynomial_degree + 1), M>=1 | |
| | | | |
| OUTPUT PARAMETERS: | | OUTPUT PARAMETERS: | |
|
| Info - error code: | | Info- same format as in LSFitLinearW() subroutine: | |
| * -4 internal SVD decomposition subroutine failed (very | | * Info>0 task is solved | |
| rare and for degenerate systems only) | | * Info<=0 an error occured: | |
| * -3 either too many constraints (M or more), | | -4 means inconvergence of internal SVD | |
| degenerate constraints (some constraints are | | -3 means inconsistent constraints | |
| repetead twice) or inconsistent constraints were | | P - interpolant in barycentric form. | |
| specified. | | Rep - report, same format as in LSFitLinearW() subroutine. | |
| * 1 task is solved | | Following fields are set: | |
| C - decomposition coefficients, array[0..M-1] | | * RMSError rms error on the (X,Y). | |
| Rep - fitting report. Following fields are set: | | * AvgError average error on the (X,Y). | |
| * RMSError rms error on the (X,Y). | | * AvgRelError average relative error on the non-zero Y | |
| * AvgError average error on the (X,Y). | | * MaxError maximum error | |
| * AvgRelError average relative error on the non-zero | | NON-WEIGHTED ERRORS ARE CALCULATED | |
| Y | | | |
| * MaxError maximum error | | | |
| NON-WEIGHTED ERRORS ARE CALCULATED | | | |
| | | | |
| IMPORTANT: | | IMPORTANT: | |
| this subroitine doesn't calculate task's condition number for K<>0. | | this subroitine doesn't calculate task's condition number for K<>0. | |
| | | | |
|
| SEE ALSO | | NOTES: | |
| LSFitLinear | | you can convert P from barycentric form to the power or Chebyshev | |
| LSFitLinearC | | basis with PolynomialBar2Pow() or PolynomialBar2Cheb() functions from | |
| LSFitLinearWC | | POLINT subpackage. | |
| | | | |
| -- ALGLIB -- | | | |
| Copyright 07.09.2009 by Bochkanov Sergey | | | |
| *************************************************************************/ | | | |
| void lsfitlinearwc(const real_1d_array &y, const real_1d_array &w, const re | | | |
| al_2d_array &fmatrix, const real_2d_array &cmatrix, const ae_int_t n, const | | | |
| ae_int_t m, const ae_int_t k, ae_int_t &info, real_1d_array &c, lsfitrepor | | | |
| t &rep); | | | |
| void lsfitlinearwc(const real_1d_array &y, const real_1d_array &w, const re | | | |
| al_2d_array &fmatrix, const real_2d_array &cmatrix, ae_int_t &info, real_1d | | | |
| _array &c, lsfitreport &rep); | | | |
| | | | |
| /************************************************************************* | | | |
| Linear least squares fitting, without weights. | | | |
| | | | |
| See LSFitLinearW for more information. | | | |
| | | | |
|
| -- ALGLIB -- | | SETTING CONSTRAINTS - DANGERS AND OPPORTUNITIES: | |
| Copyright 17.08.2009 by Bochkanov Sergey | | | |
| *************************************************************************/ | | | |
| void lsfitlinear(const real_1d_array &y, const real_2d_array &fmatrix, cons | | | |
| t ae_int_t n, const ae_int_t m, ae_int_t &info, real_1d_array &c, lsfitrepo | | | |
| rt &rep); | | | |
| void lsfitlinear(const real_1d_array &y, const real_2d_array &fmatrix, ae_i | | | |
| nt_t &info, real_1d_array &c, lsfitreport &rep); | | | |
| | | | |
|
| /************************************************************************* | | Setting constraints can lead to undesired results, like ill-conditioned | |
| Constained linear least squares fitting, without weights. | | behavior, or inconsistency being detected. From the other side, it allows | |
| | | us to improve quality of the fit. Here we summarize our experience with | |
| | | constrained regression splines: | |
| | | * even simple constraints can be inconsistent, see Wikipedia article on | |
| | | this subject: http://en.wikipedia.org/wiki/Birkhoff_interpolation | |
| | | * the greater is M (given fixed constraints), the more chances that | |
| | | constraints will be consistent | |
| | | * in the general case, consistency of constraints is NOT GUARANTEED. | |
| | | * in the one special cases, however, we can guarantee consistency. This | |
| | | case is: M>1 and constraints on the function values (NOT DERIVATIVES) | |
| | | | |
|
| See LSFitLinearWC() for more information. | | Our final recommendation is to use constraints WHEN AND ONLY when you | |
| | | can't solve your task without them. Anything beyond special cases given | |
| | | above is not guaranteed and may result in inconsistency. | |
| | | | |
|
| -- ALGLIB -- | | -- ALGLIB PROJECT -- | |
| Copyright 07.09.2009 by Bochkanov Sergey | | Copyright 10.12.2009 by Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
|
| void lsfitlinearc(const real_1d_array &y, const real_2d_array &fmatrix, con | | void polynomialfitwc(const real_1d_array &x, const real_1d_array &y, const | |
| st real_2d_array &cmatrix, const ae_int_t n, const ae_int_t m, const ae_int | | real_1d_array &w, const ae_int_t n, const real_1d_array &xc, const real_1d_ | |
| _t k, ae_int_t &info, real_1d_array &c, lsfitreport &rep); | | array &yc, const integer_1d_array &dc, const ae_int_t k, const ae_int_t m, | |
| void lsfitlinearc(const real_1d_array &y, const real_2d_array &fmatrix, con | | ae_int_t &info, barycentricinterpolant &p, polynomialfitreport &rep); | |
| st real_2d_array &cmatrix, ae_int_t &info, real_1d_array &c, lsfitreport &r | | void polynomialfitwc(const real_1d_array &x, const real_1d_array &y, const | |
| ep); | | real_1d_array &w, const real_1d_array &xc, const real_1d_array &yc, const i | |
| | | nteger_1d_array &dc, const ae_int_t m, ae_int_t &info, barycentricinterpola | |
| | | nt &p, polynomialfitreport &rep); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
|
| Weighted nonlinear least squares fitting using gradient only. | | Weghted rational least squares fitting using Floater-Hormann rational | |
| | | functions with optimal D chosen from [0,9], with constraints and | |
| Nonlinear task min(F(c)) is solved, where | | individual weights. | |
| | | | |
| F(c) = (w[0]*(f(c,x[0])-y[0]))^2 + ... + (w[n-1]*(f(c,x[n-1])-y[n-1]))^ | | | |
| 2, | | | |
| | | | |
|
| * N is a number of points, | | Equidistant grid with M node on [min(x),max(x)] is used to build basis | |
| * M is a dimension of a space points belong to, | | functions. Different values of D are tried, optimal D (least WEIGHTED root | |
| * K is a dimension of a space of parameters being fitted, | | mean square error) is chosen. Task is linear, so linear least squares | |
| * w is an N-dimensional vector of weight coefficients, | | solver is used. Complexity of this computational scheme is O(N*M^2) | |
| * x is a set of N points, each of them is an M-dimensional vector, | | (mostly dominated by the least squares solver). | |
| * c is a K-dimensional vector of parameters being fitted | | | |
| | | | |
|
| This subroutine uses only f(c,x[i]) and its gradient. | | SEE ALSO | |
| | | * BarycentricFitFloaterHormann(), "lightweight" fitting without invididual | |
| | | weights and constraints. | |
| | | | |
| INPUT PARAMETERS: | | INPUT PARAMETERS: | |
|
| X - array[0..N-1,0..M-1], points (one row = one point) | | X - points, array[0..N-1]. | |
| Y - array[0..N-1], function values. | | Y - function values, array[0..N-1]. | |
| W - weights, array[0..N-1] | | W - weights, array[0..N-1] | |
| C - array[0..K-1], initial approximation to the solution, | | Each summand in square sum of approximation deviations from | |
| N - number of points, N>1 | | given values is multiplied by the square of corresponding | |
| M - dimension of space | | weight. Fill it by 1's if you don't want to solve weighted | |
| | | task. | |
| | | N - number of points, N>0. | |
| | | XC - points where function values/derivatives are constrained, | |
| | | array[0..K-1]. | |
| | | YC - values of constraints, array[0..K-1] | |
| | | DC - array[0..K-1], types of constraints: | |
| | | * DC[i]=0 means that S(XC[i])=YC[i] | |
| | | * DC[i]=1 means that S'(XC[i])=YC[i] | |
| | | SEE BELOW FOR IMPORTANT INFORMATION ON CONSTRAINTS | |
| | | K - number of constraints, 0<=K<M. | |
| | | K=0 means no constraints (XC/YC/DC are not used in such cases) | |
| | | M - number of basis functions ( = number_of_nodes), M>=2. | |
| | | | |
| | | OUTPUT PARAMETERS: | |
| | | Info- same format as in LSFitLinearWC() subroutine. | |
| | | * Info>0 task is solved | |
| | | * Info<=0 an error occured: | |
| | | -4 means inconvergence of internal SVD | |
| | | -3 means inconsistent constraints | |
| | | -1 means another errors in parameters passed | |
| | | (N<=0, for example) | |
| | | B - barycentric interpolant. | |
| | | Rep - report, same format as in LSFitLinearWC() subroutine. | |
| | | Following fields are set: | |
| | | * DBest best value of the D parameter | |
| | | * RMSError rms error on the (X,Y). | |
| | | * AvgError average error on the (X,Y). | |
| | | * AvgRelError average relative error on the non-zero Y | |
| | | * MaxError maximum error | |
| | | NON-WEIGHTED ERRORS ARE CALCULATED | |
| | | | |
| | | IMPORTANT: | |
| | | this subroutine doesn't calculate task's condition number for K<>0. | |
| | | | |
| | | SETTING CONSTRAINTS - DANGERS AND OPPORTUNITIES: | |
| | | | |
| | | Setting constraints can lead to undesired results, like ill-conditioned | |
| | | behavior, or inconsistency being detected. From the other side, it allows | |
| | | us to improve quality of the fit. Here we summarize our experience with | |
| | | constrained barycentric interpolants: | |
| | | * excessive constraints can be inconsistent. Floater-Hormann basis | |
| | | functions aren't as flexible as splines (although they are very smooth). | |
| | | * the more evenly constraints are spread across [min(x),max(x)], the more | |
| | | chances that they will be consistent | |
| | | * the greater is M (given fixed constraints), the more chances that | |
| | | constraints will be consistent | |
| | | * in the general case, consistency of constraints IS NOT GUARANTEED. | |
| | | * in the several special cases, however, we CAN guarantee consistency. | |
| | | * one of this cases is constraints on the function VALUES at the interval | |
| | | boundaries. Note that consustency of the constraints on the function | |
| | | DERIVATIVES is NOT guaranteed (you can use in such cases cubic splines | |
| | | which are more flexible). | |
| | | * another special case is ONE constraint on the function value (OR, but | |
| | | not AND, derivative) anywhere in the interval | |
| | | | |
| | | Our final recommendation is to use constraints WHEN AND ONLY WHEN you | |
| | | can't solve your task without them. Anything beyond special cases given | |
| | | above is not guaranteed and may result in inconsistency. | |
| | | | |
| | | -- ALGLIB PROJECT -- | |
| | | Copyright 18.08.2009 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void barycentricfitfloaterhormannwc(const real_1d_array &x, const real_1d_a | |
| | | rray &y, const real_1d_array &w, const ae_int_t n, const real_1d_array &xc, | |
| | | const real_1d_array &yc, const integer_1d_array &dc, const ae_int_t k, con | |
| | | st ae_int_t m, ae_int_t &info, barycentricinterpolant &b, barycentricfitrep | |
| | | ort &rep); | |
| | | | |
| | | /************************************************************************* | |
| | | Rational least squares fitting using Floater-Hormann rational functions | |
| | | with optimal D chosen from [0,9]. | |
| | | | |
| | | Equidistant grid with M node on [min(x),max(x)] is used to build basis | |
| | | functions. Different values of D are tried, optimal D (least root mean | |
| | | square error) is chosen. Task is linear, so linear least squares solver | |
| | | is used. Complexity of this computational scheme is O(N*M^2) (mostly | |
| | | dominated by the least squares solver). | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | X - points, array[0..N-1]. | |
| | | Y - function values, array[0..N-1]. | |
| | | N - number of points, N>0. | |
| | | M - number of basis functions ( = number_of_nodes), M>=2. | |
| | | | |
| | | OUTPUT PARAMETERS: | |
| | | Info- same format as in LSFitLinearWC() subroutine. | |
| | | * Info>0 task is solved | |
| | | * Info<=0 an error occured: | |
| | | -4 means inconvergence of internal SVD | |
| | | -3 means inconsistent constraints | |
| | | B - barycentric interpolant. | |
| | | Rep - report, same format as in LSFitLinearWC() subroutine. | |
| | | Following fields are set: | |
| | | * DBest best value of the D parameter | |
| | | * RMSError rms error on the (X,Y). | |
| | | * AvgError average error on the (X,Y). | |
| | | * AvgRelError average relative error on the non-zero Y | |
| | | * MaxError maximum error | |
| | | NON-WEIGHTED ERRORS ARE CALCULATED | |
| | | | |
| | | -- ALGLIB PROJECT -- | |
| | | Copyright 18.08.2009 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void barycentricfitfloaterhormann(const real_1d_array &x, const real_1d_arr | |
| | | ay &y, const ae_int_t n, const ae_int_t m, ae_int_t &info, barycentricinter | |
| | | polant &b, barycentricfitreport &rep); | |
| | | | |
| | | /************************************************************************* | |
| | | Rational least squares fitting using Floater-Hormann rational functions | |
| | | with optimal D chosen from [0,9]. | |
| | | | |
| | | Equidistant grid with M node on [min(x),max(x)] is used to build basis | |
| | | functions. Different values of D are tried, optimal D (least root mean | |
| | | square error) is chosen. Task is linear, so linear least squares solver | |
| | | is used. Complexity of this computational scheme is O(N*M^2) (mostly | |
| | | dominated by the least squares solver). | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | X - points, array[0..N-1]. | |
| | | Y - function values, array[0..N-1]. | |
| | | N - number of points, N>0. | |
| | | M - number of basis functions ( = number_of_nodes), M>=2. | |
| | | | |
| | | OUTPUT PARAMETERS: | |
| | | Info- same format as in LSFitLinearWC() subroutine. | |
| | | * Info>0 task is solved | |
| | | * Info<=0 an error occured: | |
| | | -4 means inconvergence of internal SVD | |
| | | -3 means inconsistent constraints | |
| | | B - barycentric interpolant. | |
| | | Rep - report, same format as in LSFitLinearWC() subroutine. | |
| | | Following fields are set: | |
| | | * DBest best value of the D parameter | |
| | | * RMSError rms error on the (X,Y). | |
| | | * AvgError average error on the (X,Y). | |
| | | * AvgRelError average relative error on the non-zero Y | |
| | | * MaxError maximum error | |
| | | NON-WEIGHTED ERRORS ARE CALCULATED | |
| | | | |
| | | -- ALGLIB PROJECT -- | |
| | | Copyright 18.08.2009 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void spline1dfitpenalized(const real_1d_array &x, const real_1d_array &y, c | |
| | | onst ae_int_t n, const ae_int_t m, const double rho, ae_int_t &info, spline | |
| | | 1dinterpolant &s, spline1dfitreport &rep); | |
| | | void spline1dfitpenalized(const real_1d_array &x, const real_1d_array &y, c | |
| | | onst ae_int_t m, const double rho, ae_int_t &info, spline1dinterpolant &s, | |
| | | spline1dfitreport &rep); | |
| | | | |
| | | /************************************************************************* | |
| | | Weighted fitting by penalized cubic spline. | |
| | | | |
| | | Equidistant grid with M nodes on [min(x,xc),max(x,xc)] is used to build | |
| | | basis functions. Basis functions are cubic splines with natural boundary | |
| | | conditions. Problem is regularized by adding non-linearity penalty to the | |
| | | usual least squares penalty function: | |
| | | | |
| | | S(x) = arg min { LS + P }, where | |
| | | LS = SUM { w[i]^2*(y[i] - S(x[i]))^2 } - least squares penalty | |
| | | P = C*10^rho*integral{ S''(x)^2*dx } - non-linearity penalty | |
| | | rho - tunable constant given by user | |
| | | C - automatically determined scale parameter, | |
| | | makes penalty invariant with respect to scaling of X, Y, W. | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | X - points, array[0..N-1]. | |
| | | Y - function values, array[0..N-1]. | |
| | | W - weights, array[0..N-1] | |
| | | Each summand in square sum of approximation deviations from | |
| | | given values is multiplied by the square of corresponding | |
| | | weight. Fill it by 1's if you don't want to solve weighted | |
| | | problem. | |
| | | N - number of points (optional): | |
| | | * N>0 | |
| | | * if given, only first N elements of X/Y/W are processed | |
| | | * if not given, automatically determined from X/Y/W sizes | |
| | | M - number of basis functions ( = number_of_nodes), M>=4. | |
| | | Rho - regularization constant passed by user. It penalizes | |
| | | nonlinearity in the regression spline. It is logarithmically | |
| | | scaled, i.e. actual value of regularization constant is | |
| | | calculated as 10^Rho. It is automatically scaled so that: | |
| | | * Rho=2.0 corresponds to moderate amount of nonlinearity | |
| | | * generally, it should be somewhere in the [-8.0,+8.0] | |
| | | If you do not want to penalize nonlineary, | |
| | | pass small Rho. Values as low as -15 should work. | |
| | | | |
| | | OUTPUT PARAMETERS: | |
| | | Info- same format as in LSFitLinearWC() subroutine. | |
| | | * Info>0 task is solved | |
| | | * Info<=0 an error occured: | |
| | | -4 means inconvergence of internal SVD or | |
| | | Cholesky decomposition; problem may be | |
| | | too ill-conditioned (very rare) | |
| | | S - spline interpolant. | |
| | | Rep - Following fields are set: | |
| | | * RMSError rms error on the (X,Y). | |
| | | * AvgError average error on the (X,Y). | |
| | | * AvgRelError average relative error on the non-zero Y | |
| | | * MaxError maximum error | |
| | | NON-WEIGHTED ERRORS ARE CALCULATED | |
| | | | |
| | | IMPORTANT: | |
| | | this subroitine doesn't calculate task's condition number for K<>0. | |
| | | | |
| | | NOTE 1: additional nodes are added to the spline outside of the fitting | |
| | | interval to force linearity when x<min(x,xc) or x>max(x,xc). It is done | |
| | | for consistency - we penalize non-linearity at [min(x,xc),max(x,xc)], so | |
| | | it is natural to force linearity outside of this interval. | |
| | | | |
| | | NOTE 2: function automatically sorts points, so caller may pass unsorted | |
| | | array. | |
| | | | |
| | | -- ALGLIB PROJECT -- | |
| | | Copyright 19.10.2010 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void spline1dfitpenalizedw(const real_1d_array &x, const real_1d_array &y, | |
| | | const real_1d_array &w, const ae_int_t n, const ae_int_t m, const double rh | |
| | | o, ae_int_t &info, spline1dinterpolant &s, spline1dfitreport &rep); | |
| | | void spline1dfitpenalizedw(const real_1d_array &x, const real_1d_array &y, | |
| | | const real_1d_array &w, const ae_int_t m, const double rho, ae_int_t &info, | |
| | | spline1dinterpolant &s, spline1dfitreport &rep); | |
| | | | |
| | | /************************************************************************* | |
| | | Weighted fitting by cubic spline, with constraints on function values or | |
| | | derivatives. | |
| | | | |
| | | Equidistant grid with M-2 nodes on [min(x,xc),max(x,xc)] is used to build | |
| | | basis functions. Basis functions are cubic splines with continuous second | |
| | | derivatives and non-fixed first derivatives at interval ends. Small | |
| | | regularizing term is used when solving constrained tasks (to improve | |
| | | stability). | |
| | | | |
| | | Task is linear, so linear least squares solver is used. Complexity of this | |
| | | computational scheme is O(N*M^2), mostly dominated by least squares solver | |
| | | | |
| | | SEE ALSO | |
| | | Spline1DFitHermiteWC() - fitting by Hermite splines (more flexible, | |
| | | less smooth) | |
| | | Spline1DFitCubic() - "lightweight" fitting by cubic splines, | |
| | | without invididual weights and constraints | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | X - points, array[0..N-1]. | |
| | | Y - function values, array[0..N-1]. | |
| | | W - weights, array[0..N-1] | |
| | | Each summand in square sum of approximation deviations from | |
| | | given values is multiplied by the square of corresponding | |
| | | weight. Fill it by 1's if you don't want to solve weighted | |
| | | task. | |
| | | N - number of points (optional): | |
| | | * N>0 | |
| | | * if given, only first N elements of X/Y/W are processed | |
| | | * if not given, automatically determined from X/Y/W sizes | |
| | | XC - points where spline values/derivatives are constrained, | |
| | | array[0..K-1]. | |
| | | YC - values of constraints, array[0..K-1] | |
| | | DC - array[0..K-1], types of constraints: | |
| | | * DC[i]=0 means that S(XC[i])=YC[i] | |
| | | * DC[i]=1 means that S'(XC[i])=YC[i] | |
| | | SEE BELOW FOR IMPORTANT INFORMATION ON CONSTRAINTS | |
| | | K - number of constraints (optional): | |
| | | * 0<=K<M. | |
| | | * K=0 means no constraints (XC/YC/DC are not used) | |
| | | * if given, only first K elements of XC/YC/DC are used | |
| | | * if not given, automatically determined from XC/YC/DC | |
| | | M - number of basis functions ( = number_of_nodes+2), M>=4. | |
| | | | |
| | | OUTPUT PARAMETERS: | |
| | | Info- same format as in LSFitLinearWC() subroutine. | |
| | | * Info>0 task is solved | |
| | | * Info<=0 an error occured: | |
| | | -4 means inconvergence of internal SVD | |
| | | -3 means inconsistent constraints | |
| | | S - spline interpolant. | |
| | | Rep - report, same format as in LSFitLinearWC() subroutine. | |
| | | Following fields are set: | |
| | | * RMSError rms error on the (X,Y). | |
| | | * AvgError average error on the (X,Y). | |
| | | * AvgRelError average relative error on the non-zero Y | |
| | | * MaxError maximum error | |
| | | NON-WEIGHTED ERRORS ARE CALCULATED | |
| | | | |
| | | IMPORTANT: | |
| | | this subroitine doesn't calculate task's condition number for K<>0. | |
| | | | |
| | | ORDER OF POINTS | |
| | | | |
| | | Subroutine automatically sorts points, so caller may pass unsorted array. | |
| | | | |
| | | SETTING CONSTRAINTS - DANGERS AND OPPORTUNITIES: | |
| | | | |
| | | Setting constraints can lead to undesired results, like ill-conditioned | |
| | | behavior, or inconsistency being detected. From the other side, it allows | |
| | | us to improve quality of the fit. Here we summarize our experience with | |
| | | constrained regression splines: | |
| | | * excessive constraints can be inconsistent. Splines are piecewise cubic | |
| | | functions, and it is easy to create an example, where large number of | |
| | | constraints concentrated in small area will result in inconsistency. | |
| | | Just because spline is not flexible enough to satisfy all of them. And | |
| | | same constraints spread across the [min(x),max(x)] will be perfectly | |
| | | consistent. | |
| | | * the more evenly constraints are spread across [min(x),max(x)], the more | |
| | | chances that they will be consistent | |
| | | * the greater is M (given fixed constraints), the more chances that | |
| | | constraints will be consistent | |
| | | * in the general case, consistency of constraints IS NOT GUARANTEED. | |
| | | * in the several special cases, however, we CAN guarantee consistency. | |
| | | * one of this cases is constraints on the function values AND/OR its | |
| | | derivatives at the interval boundaries. | |
| | | * another special case is ONE constraint on the function value (OR, but | |
| | | not AND, derivative) anywhere in the interval | |
| | | | |
| | | Our final recommendation is to use constraints WHEN AND ONLY WHEN you | |
| | | can't solve your task without them. Anything beyond special cases given | |
| | | above is not guaranteed and may result in inconsistency. | |
| | | | |
| | | -- ALGLIB PROJECT -- | |
| | | Copyright 18.08.2009 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void spline1dfitcubicwc(const real_1d_array &x, const real_1d_array &y, con | |
| | | st real_1d_array &w, const ae_int_t n, const real_1d_array &xc, const real_ | |
| | | 1d_array &yc, const integer_1d_array &dc, const ae_int_t k, const ae_int_t | |
| | | m, ae_int_t &info, spline1dinterpolant &s, spline1dfitreport &rep); | |
| | | void spline1dfitcubicwc(const real_1d_array &x, const real_1d_array &y, con | |
| | | st real_1d_array &w, const real_1d_array &xc, const real_1d_array &yc, cons | |
| | | t integer_1d_array &dc, const ae_int_t m, ae_int_t &info, spline1dinterpola | |
| | | nt &s, spline1dfitreport &rep); | |
| | | | |
| | | /************************************************************************* | |
| | | Weighted fitting by Hermite spline, with constraints on function values | |
| | | or first derivatives. | |
| | | | |
| | | Equidistant grid with M nodes on [min(x,xc),max(x,xc)] is used to build | |
| | | basis functions. Basis functions are Hermite splines. Small regularizing | |
| | | term is used when solving constrained tasks (to improve stability). | |
| | | | |
| | | Task is linear, so linear least squares solver is used. Complexity of this | |
| | | computational scheme is O(N*M^2), mostly dominated by least squares solver | |
| | | | |
| | | SEE ALSO | |
| | | Spline1DFitCubicWC() - fitting by Cubic splines (less flexible, | |
| | | more smooth) | |
| | | Spline1DFitHermite() - "lightweight" Hermite fitting, without | |
| | | invididual weights and constraints | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | X - points, array[0..N-1]. | |
| | | Y - function values, array[0..N-1]. | |
| | | W - weights, array[0..N-1] | |
| | | Each summand in square sum of approximation deviations from | |
| | | given values is multiplied by the square of corresponding | |
| | | weight. Fill it by 1's if you don't want to solve weighted | |
| | | task. | |
| | | N - number of points (optional): | |
| | | * N>0 | |
| | | * if given, only first N elements of X/Y/W are processed | |
| | | * if not given, automatically determined from X/Y/W sizes | |
| | | XC - points where spline values/derivatives are constrained, | |
| | | array[0..K-1]. | |
| | | YC - values of constraints, array[0..K-1] | |
| | | DC - array[0..K-1], types of constraints: | |
| | | * DC[i]=0 means that S(XC[i])=YC[i] | |
| | | * DC[i]=1 means that S'(XC[i])=YC[i] | |
| | | SEE BELOW FOR IMPORTANT INFORMATION ON CONSTRAINTS | |
| | | K - number of constraints (optional): | |
| | | * 0<=K<M. | |
| | | * K=0 means no constraints (XC/YC/DC are not used) | |
| | | * if given, only first K elements of XC/YC/DC are used | |
| | | * if not given, automatically determined from XC/YC/DC | |
| | | M - number of basis functions (= 2 * number of nodes), | |
| | | M>=4, | |
| | | M IS EVEN! | |
| | | | |
| | | OUTPUT PARAMETERS: | |
| | | Info- same format as in LSFitLinearW() subroutine: | |
| | | * Info>0 task is solved | |
| | | * Info<=0 an error occured: | |
| | | -4 means inconvergence of internal SVD | |
| | | -3 means inconsistent constraints | |
| | | -2 means odd M was passed (which is not supported) | |
| | | -1 means another errors in parameters passed | |
| | | (N<=0, for example) | |
| | | S - spline interpolant. | |
| | | Rep - report, same format as in LSFitLinearW() subroutine. | |
| | | Following fields are set: | |
| | | * RMSError rms error on the (X,Y). | |
| | | * AvgError average error on the (X,Y). | |
| | | * AvgRelError average relative error on the non-zero Y | |
| | | * MaxError maximum error | |
| | | NON-WEIGHTED ERRORS ARE CALCULATED | |
| | | | |
| | | IMPORTANT: | |
| | | this subroitine doesn't calculate task's condition number for K<>0. | |
| | | | |
| | | IMPORTANT: | |
| | | this subroitine supports only even M's | |
| | | | |
| | | ORDER OF POINTS | |
| | | | |
| | | Subroutine automatically sorts points, so caller may pass unsorted array. | |
| | | | |
| | | SETTING CONSTRAINTS - DANGERS AND OPPORTUNITIES: | |
| | | | |
| | | Setting constraints can lead to undesired results, like ill-conditioned | |
| | | behavior, or inconsistency being detected. From the other side, it allows | |
| | | us to improve quality of the fit. Here we summarize our experience with | |
| | | constrained regression splines: | |
| | | * excessive constraints can be inconsistent. Splines are piecewise cubic | |
| | | functions, and it is easy to create an example, where large number of | |
| | | constraints concentrated in small area will result in inconsistency. | |
| | | Just because spline is not flexible enough to satisfy all of them. And | |
| | | same constraints spread across the [min(x),max(x)] will be perfectly | |
| | | consistent. | |
| | | * the more evenly constraints are spread across [min(x),max(x)], the more | |
| | | chances that they will be consistent | |
| | | * the greater is M (given fixed constraints), the more chances that | |
| | | constraints will be consistent | |
| | | * in the general case, consistency of constraints is NOT GUARANTEED. | |
| | | * in the several special cases, however, we can guarantee consistency. | |
| | | * one of this cases is M>=4 and constraints on the function value | |
| | | (AND/OR its derivative) at the interval boundaries. | |
| | | * another special case is M>=4 and ONE constraint on the function value | |
| | | (OR, BUT NOT AND, derivative) anywhere in [min(x),max(x)] | |
| | | | |
| | | Our final recommendation is to use constraints WHEN AND ONLY when you | |
| | | can't solve your task without them. Anything beyond special cases given | |
| | | above is not guaranteed and may result in inconsistency. | |
| | | | |
| | | -- ALGLIB PROJECT -- | |
| | | Copyright 18.08.2009 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void spline1dfithermitewc(const real_1d_array &x, const real_1d_array &y, c | |
| | | onst real_1d_array &w, const ae_int_t n, const real_1d_array &xc, const rea | |
| | | l_1d_array &yc, const integer_1d_array &dc, const ae_int_t k, const ae_int_ | |
| | | t m, ae_int_t &info, spline1dinterpolant &s, spline1dfitreport &rep); | |
| | | void spline1dfithermitewc(const real_1d_array &x, const real_1d_array &y, c | |
| | | onst real_1d_array &w, const real_1d_array &xc, const real_1d_array &yc, co | |
| | | nst integer_1d_array &dc, const ae_int_t m, ae_int_t &info, spline1dinterpo | |
| | | lant &s, spline1dfitreport &rep); | |
| | | | |
| | | /************************************************************************* | |
| | | Least squares fitting by cubic spline. | |
| | | | |
| | | This subroutine is "lightweight" alternative for more complex and feature- | |
| | | rich Spline1DFitCubicWC(). See Spline1DFitCubicWC() for more information | |
| | | about subroutine parameters (we don't duplicate it here because of length) | |
| | | | |
| | | -- ALGLIB PROJECT -- | |
| | | Copyright 18.08.2009 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void spline1dfitcubic(const real_1d_array &x, const real_1d_array &y, const | |
| | | ae_int_t n, const ae_int_t m, ae_int_t &info, spline1dinterpolant &s, spli | |
| | | ne1dfitreport &rep); | |
| | | void spline1dfitcubic(const real_1d_array &x, const real_1d_array &y, const | |
| | | ae_int_t m, ae_int_t &info, spline1dinterpolant &s, spline1dfitreport &rep | |
| | | ); | |
| | | | |
| | | /************************************************************************* | |
| | | Least squares fitting by Hermite spline. | |
| | | | |
| | | This subroutine is "lightweight" alternative for more complex and feature- | |
| | | rich Spline1DFitHermiteWC(). See Spline1DFitHermiteWC() description for | |
| | | more information about subroutine parameters (we don't duplicate it here | |
| | | because of length). | |
| | | | |
| | | -- ALGLIB PROJECT -- | |
| | | Copyright 18.08.2009 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void spline1dfithermite(const real_1d_array &x, const real_1d_array &y, con | |
| | | st ae_int_t n, const ae_int_t m, ae_int_t &info, spline1dinterpolant &s, sp | |
| | | line1dfitreport &rep); | |
| | | void spline1dfithermite(const real_1d_array &x, const real_1d_array &y, con | |
| | | st ae_int_t m, ae_int_t &info, spline1dinterpolant &s, spline1dfitreport &r | |
| | | ep); | |
| | | | |
| | | /************************************************************************* | |
| | | Weighted linear least squares fitting. | |
| | | | |
| | | QR decomposition is used to reduce task to MxM, then triangular solver or | |
| | | SVD-based solver is used depending on condition number of the system. It | |
| | | allows to maximize speed and retain decent accuracy. | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | Y - array[0..N-1] Function values in N points. | |
| | | W - array[0..N-1] Weights corresponding to function values. | |
| | | Each summand in square sum of approximation deviations | |
| | | from given values is multiplied by the square of | |
| | | corresponding weight. | |
| | | FMatrix - a table of basis functions values, array[0..N-1, 0..M-1]. | |
| | | FMatrix[I, J] - value of J-th basis function in I-th point. | |
| | | N - number of points used. N>=1. | |
| | | M - number of basis functions, M>=1. | |
| | | | |
| | | OUTPUT PARAMETERS: | |
| | | Info - error code: | |
| | | * -4 internal SVD decomposition subroutine failed (very | |
| | | rare and for degenerate systems only) | |
| | | * -1 incorrect N/M were specified | |
| | | * 1 task is solved | |
| | | C - decomposition coefficients, array[0..M-1] | |
| | | Rep - fitting report. Following fields are set: | |
| | | * Rep.TaskRCond reciprocal of condition number | |
| | | * RMSError rms error on the (X,Y). | |
| | | * AvgError average error on the (X,Y). | |
| | | * AvgRelError average relative error on the non-zero | |
| | | Y | |
| | | * MaxError maximum error | |
| | | NON-WEIGHTED ERRORS ARE CALCULATED | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 17.08.2009 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void lsfitlinearw(const real_1d_array &y, const real_1d_array &w, const rea | |
| | | l_2d_array &fmatrix, const ae_int_t n, const ae_int_t m, ae_int_t &info, re | |
| | | al_1d_array &c, lsfitreport &rep); | |
| | | void lsfitlinearw(const real_1d_array &y, const real_1d_array &w, const rea | |
| | | l_2d_array &fmatrix, ae_int_t &info, real_1d_array &c, lsfitreport &rep); | |
| | | | |
| | | /************************************************************************* | |
| | | Weighted constained linear least squares fitting. | |
| | | | |
| | | This is variation of LSFitLinearW(), which searchs for min|A*x=b| given | |
| | | that K additional constaints C*x=bc are satisfied. It reduces original | |
| | | task to modified one: min|B*y-d| WITHOUT constraints, then LSFitLinearW() | |
| | | is called. | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | Y - array[0..N-1] Function values in N points. | |
| | | W - array[0..N-1] Weights corresponding to function values. | |
| | | Each summand in square sum of approximation deviations | |
| | | from given values is multiplied by the square of | |
| | | corresponding weight. | |
| | | FMatrix - a table of basis functions values, array[0..N-1, 0..M-1]. | |
| | | FMatrix[I,J] - value of J-th basis function in I-th point. | |
| | | CMatrix - a table of constaints, array[0..K-1,0..M]. | |
| | | I-th row of CMatrix corresponds to I-th linear constraint: | |
| | | CMatrix[I,0]*C[0] + ... + CMatrix[I,M-1]*C[M-1] = CMatrix[I | |
| | | ,M] | |
| | | N - number of points used. N>=1. | |
| | | M - number of basis functions, M>=1. | |
| | | K - number of constraints, 0 <= K < M | |
| | | K=0 corresponds to absence of constraints. | |
| | | | |
| | | OUTPUT PARAMETERS: | |
| | | Info - error code: | |
| | | * -4 internal SVD decomposition subroutine failed (very | |
| | | rare and for degenerate systems only) | |
| | | * -3 either too many constraints (M or more), | |
| | | degenerate constraints (some constraints are | |
| | | repetead twice) or inconsistent constraints were | |
| | | specified. | |
| | | * 1 task is solved | |
| | | C - decomposition coefficients, array[0..M-1] | |
| | | Rep - fitting report. Following fields are set: | |
| | | * RMSError rms error on the (X,Y). | |
| | | * AvgError average error on the (X,Y). | |
| | | * AvgRelError average relative error on the non-zero | |
| | | Y | |
| | | * MaxError maximum error | |
| | | NON-WEIGHTED ERRORS ARE CALCULATED | |
| | | | |
| | | IMPORTANT: | |
| | | this subroitine doesn't calculate task's condition number for K<>0. | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 07.09.2009 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void lsfitlinearwc(const real_1d_array &y, const real_1d_array &w, const re | |
| | | al_2d_array &fmatrix, const real_2d_array &cmatrix, const ae_int_t n, const | |
| | | ae_int_t m, const ae_int_t k, ae_int_t &info, real_1d_array &c, lsfitrepor | |
| | | t &rep); | |
| | | void lsfitlinearwc(const real_1d_array &y, const real_1d_array &w, const re | |
| | | al_2d_array &fmatrix, const real_2d_array &cmatrix, ae_int_t &info, real_1d | |
| | | _array &c, lsfitreport &rep); | |
| | | | |
| | | /************************************************************************* | |
| | | Linear least squares fitting. | |
| | | | |
| | | QR decomposition is used to reduce task to MxM, then triangular solver or | |
| | | SVD-based solver is used depending on condition number of the system. It | |
| | | allows to maximize speed and retain decent accuracy. | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | Y - array[0..N-1] Function values in N points. | |
| | | FMatrix - a table of basis functions values, array[0..N-1, 0..M-1]. | |
| | | FMatrix[I, J] - value of J-th basis function in I-th point. | |
| | | N - number of points used. N>=1. | |
| | | M - number of basis functions, M>=1. | |
| | | | |
| | | OUTPUT PARAMETERS: | |
| | | Info - error code: | |
| | | * -4 internal SVD decomposition subroutine failed (very | |
| | | rare and for degenerate systems only) | |
| | | * 1 task is solved | |
| | | C - decomposition coefficients, array[0..M-1] | |
| | | Rep - fitting report. Following fields are set: | |
| | | * Rep.TaskRCond reciprocal of condition number | |
| | | * RMSError rms error on the (X,Y). | |
| | | * AvgError average error on the (X,Y). | |
| | | * AvgRelError average relative error on the non-zero | |
| | | Y | |
| | | * MaxError maximum error | |
| | | NON-WEIGHTED ERRORS ARE CALCULATED | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 17.08.2009 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void lsfitlinear(const real_1d_array &y, const real_2d_array &fmatrix, cons | |
| | | t ae_int_t n, const ae_int_t m, ae_int_t &info, real_1d_array &c, lsfitrepo | |
| | | rt &rep); | |
| | | void lsfitlinear(const real_1d_array &y, const real_2d_array &fmatrix, ae_i | |
| | | nt_t &info, real_1d_array &c, lsfitreport &rep); | |
| | | | |
| | | /************************************************************************* | |
| | | Constained linear least squares fitting. | |
| | | | |
| | | This is variation of LSFitLinear(), which searchs for min|A*x=b| given | |
| | | that K additional constaints C*x=bc are satisfied. It reduces original | |
| | | task to modified one: min|B*y-d| WITHOUT constraints, then LSFitLinear() | |
| | | is called. | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | Y - array[0..N-1] Function values in N points. | |
| | | FMatrix - a table of basis functions values, array[0..N-1, 0..M-1]. | |
| | | FMatrix[I,J] - value of J-th basis function in I-th point. | |
| | | CMatrix - a table of constaints, array[0..K-1,0..M]. | |
| | | I-th row of CMatrix corresponds to I-th linear constraint: | |
| | | CMatrix[I,0]*C[0] + ... + CMatrix[I,M-1]*C[M-1] = CMatrix[I | |
| | | ,M] | |
| | | N - number of points used. N>=1. | |
| | | M - number of basis functions, M>=1. | |
| | | K - number of constraints, 0 <= K < M | |
| | | K=0 corresponds to absence of constraints. | |
| | | | |
| | | OUTPUT PARAMETERS: | |
| | | Info - error code: | |
| | | * -4 internal SVD decomposition subroutine failed (very | |
| | | rare and for degenerate systems only) | |
| | | * -3 either too many constraints (M or more), | |
| | | degenerate constraints (some constraints are | |
| | | repetead twice) or inconsistent constraints were | |
| | | specified. | |
| | | * 1 task is solved | |
| | | C - decomposition coefficients, array[0..M-1] | |
| | | Rep - fitting report. Following fields are set: | |
| | | * RMSError rms error on the (X,Y). | |
| | | * AvgError average error on the (X,Y). | |
| | | * AvgRelError average relative error on the non-zero | |
| | | Y | |
| | | * MaxError maximum error | |
| | | NON-WEIGHTED ERRORS ARE CALCULATED | |
| | | | |
| | | IMPORTANT: | |
| | | this subroitine doesn't calculate task's condition number for K<>0. | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 07.09.2009 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void lsfitlinearc(const real_1d_array &y, const real_2d_array &fmatrix, con | |
| | | st real_2d_array &cmatrix, const ae_int_t n, const ae_int_t m, const ae_int | |
| | | _t k, ae_int_t &info, real_1d_array &c, lsfitreport &rep); | |
| | | void lsfitlinearc(const real_1d_array &y, const real_2d_array &fmatrix, con | |
| | | st real_2d_array &cmatrix, ae_int_t &info, real_1d_array &c, lsfitreport &r | |
| | | ep); | |
| | | | |
| | | /************************************************************************* | |
| | | Weighted nonlinear least squares fitting using function values only. | |
| | | | |
| | | Combination of numerical differentiation and secant updates is used to | |
| | | obtain function Jacobian. | |
| | | | |
| | | Nonlinear task min(F(c)) is solved, where | |
| | | | |
| | | F(c) = (w[0]*(f(c,x[0])-y[0]))^2 + ... + (w[n-1]*(f(c,x[n-1])-y[n-1]))^ | |
| | | 2, | |
| | | | |
| | | * N is a number of points, | |
| | | * M is a dimension of a space points belong to, | |
| | | * K is a dimension of a space of parameters being fitted, | |
| | | * w is an N-dimensional vector of weight coefficients, | |
| | | * x is a set of N points, each of them is an M-dimensional vector, | |
| | | * c is a K-dimensional vector of parameters being fitted | |
| | | | |
| | | This subroutine uses only f(c,x[i]). | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | X - array[0..N-1,0..M-1], points (one row = one point) | |
| | | Y - array[0..N-1], function values. | |
| | | W - weights, array[0..N-1] | |
| | | C - array[0..K-1], initial approximation to the solution, | |
| | | N - number of points, N>1 | |
| | | M - dimension of space | |
| | | K - number of parameters being fitted | |
| | | DiffStep- numerical differentiation step; | |
| | | should not be very small or large; | |
| | | large = loss of accuracy | |
| | | small = growth of round-off errors | |
| | | | |
| | | OUTPUT PARAMETERS: | |
| | | State - structure which stores algorithm state | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 18.10.2008 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void lsfitcreatewf(const real_2d_array &x, const real_1d_array &y, const re | |
| | | al_1d_array &w, const real_1d_array &c, const ae_int_t n, const ae_int_t m, | |
| | | const ae_int_t k, const double diffstep, lsfitstate &state); | |
| | | void lsfitcreatewf(const real_2d_array &x, const real_1d_array &y, const re | |
| | | al_1d_array &w, const real_1d_array &c, const double diffstep, lsfitstate & | |
| | | state); | |
| | | | |
| | | /************************************************************************* | |
| | | Nonlinear least squares fitting using function values only. | |
| | | | |
| | | Combination of numerical differentiation and secant updates is used to | |
| | | obtain function Jacobian. | |
| | | | |
| | | Nonlinear task min(F(c)) is solved, where | |
| | | | |
| | | F(c) = (f(c,x[0])-y[0])^2 + ... + (f(c,x[n-1])-y[n-1])^2, | |
| | | | |
| | | * N is a number of points, | |
| | | * M is a dimension of a space points belong to, | |
| | | * K is a dimension of a space of parameters being fitted, | |
| | | * w is an N-dimensional vector of weight coefficients, | |
| | | * x is a set of N points, each of them is an M-dimensional vector, | |
| | | * c is a K-dimensional vector of parameters being fitted | |
| | | | |
| | | This subroutine uses only f(c,x[i]). | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | X - array[0..N-1,0..M-1], points (one row = one point) | |
| | | Y - array[0..N-1], function values. | |
| | | C - array[0..K-1], initial approximation to the solution, | |
| | | N - number of points, N>1 | |
| | | M - dimension of space | |
| | | K - number of parameters being fitted | |
| | | DiffStep- numerical differentiation step; | |
| | | should not be very small or large; | |
| | | large = loss of accuracy | |
| | | small = growth of round-off errors | |
| | | | |
| | | OUTPUT PARAMETERS: | |
| | | State - structure which stores algorithm state | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 18.10.2008 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void lsfitcreatef(const real_2d_array &x, const real_1d_array &y, const rea | |
| | | l_1d_array &c, const ae_int_t n, const ae_int_t m, const ae_int_t k, const | |
| | | double diffstep, lsfitstate &state); | |
| | | void lsfitcreatef(const real_2d_array &x, const real_1d_array &y, const rea | |
| | | l_1d_array &c, const double diffstep, lsfitstate &state); | |
| | | | |
| | | /************************************************************************* | |
| | | Weighted nonlinear least squares fitting using gradient only. | |
| | | | |
| | | Nonlinear task min(F(c)) is solved, where | |
| | | | |
| | | F(c) = (w[0]*(f(c,x[0])-y[0]))^2 + ... + (w[n-1]*(f(c,x[n-1])-y[n-1]))^ | |
| | | 2, | |
| | | | |
| | | * N is a number of points, | |
| | | * M is a dimension of a space points belong to, | |
| | | * K is a dimension of a space of parameters being fitted, | |
| | | * w is an N-dimensional vector of weight coefficients, | |
| | | * x is a set of N points, each of them is an M-dimensional vector, | |
| | | * c is a K-dimensional vector of parameters being fitted | |
| | | | |
| | | This subroutine uses only f(c,x[i]) and its gradient. | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | X - array[0..N-1,0..M-1], points (one row = one point) | |
| | | Y - array[0..N-1], function values. | |
| | | W - weights, array[0..N-1] | |
| | | C - array[0..K-1], initial approximation to the solution, | |
| | | N - number of points, N>1 | |
| | | M - dimension of space | |
| K - number of parameters being fitted | | K - number of parameters being fitted | |
| CheapFG - boolean flag, which is: | | CheapFG - boolean flag, which is: | |
| * True if both function and gradient calculation complexit
y | | * True if both function and gradient calculation complexit
y | |
| are less than O(M^2). An improved algorithm can | | are less than O(M^2). An improved algorithm can | |
| be used which corresponds to FGJ scheme from | | be used which corresponds to FGJ scheme from | |
| MINLM unit. | | MINLM unit. | |
| * False otherwise. | | * False otherwise. | |
| Standard Jacibian-bases Levenberg-Marquardt algo | | Standard Jacibian-bases Levenberg-Marquardt algo | |
| will be used (FJ scheme). | | will be used (FJ scheme). | |
| | | | |
| | | | |
| skipping to change at line 864 | | skipping to change at line 1586 | |
|---|
| value func and gradient grad at given point x | | value func and gradient grad at given point x | |
| hess - callback which calculates function (or merit function) | | hess - callback which calculates function (or merit function) | |
| value func, gradient grad and Hessian hess at given point x | | value func, gradient grad and Hessian hess at given point x | |
| rep - optional callback which is called after each iteration | | rep - optional callback which is called after each iteration | |
| can be NULL | | can be NULL | |
| ptr - optional pointer which is passed to func/grad/hess/jac/rep | | ptr - optional pointer which is passed to func/grad/hess/jac/rep | |
| can be NULL | | can be NULL | |
| | | | |
| NOTES: | | NOTES: | |
| | | | |
|
| 1. this algorithm is somewhat unusual because it works with parameterized | | 1. this algorithm is somewhat unusual because it works with parameterized | |
| function f(C,X), where X is a function argument (we have many points | | function f(C,X), where X is a function argument (we have many points | |
| which are characterized by different argument values), and C is a | | which are characterized by different argument values), and C is a | |
| parameter to fit. | | parameter to fit. | |
| | | | |
| For example, if we want to do linear fit by f(c0,c1,x) = c0*x+c1, then | | | |
| x will be argument, and {c0,c1} will be parameters. | | | |
| | | | |
| It is important to understand that this algorithm finds minimum in the | | | |
| space of function PARAMETERS (not arguments), so it needs derivatives | | | |
| of f() with respect to C, not X. | | | |
| | | | |
| In the example above it will need f=c0*x+c1 and {df/dc0,df/dc1} = {x,1} | | | |
| instead of {df/dx} = {c0}. | | | |
| | | | |
| 2. Callback functions accept C as the first parameter, and X as the second | | | |
| | | | |
| 3. If state was created with LSFitCreateFG(), algorithm needs just | | | |
| function and its gradient, but if state was created with | | | |
| LSFitCreateFGH(), algorithm will need function, gradient and Hessian. | | | |
| | | | |
| According to the said above, there ase several versions of this | | | |
| function, which accept different sets of callbacks. | | | |
| | | | |
| This flexibility opens way to subtle errors - you may create state with | | | |
| LSFitCreateFGH() (optimization using Hessian), but call function which | | | |
| does not accept Hessian. So when algorithm will request Hessian, there | | | |
| will be no callback to call. In this case exception will be thrown. | | | |
| | | | |
| Be careful to avoid such errors because there is no way to find them at | | | |
| compile time - you can see them at runtime only. | | | |
| | | | |
| -- ALGLIB -- | | | |
| Copyright 17.08.2009 by Bochkanov Sergey | | | |
| | | | |
| *************************************************************************/ | | | |
| void lsfitfit(lsfitstate &state, | | | |
| void (*func)(const real_1d_array &c, const real_1d_array &x, double &fu | | | |
| nc, void *ptr), | | | |
| void (*grad)(const real_1d_array &c, const real_1d_array &x, double &fu | | | |
| nc, real_1d_array &grad, void *ptr), | | | |
| void (*rep)(const real_1d_array &c, double func, void *ptr) = NULL, | | | |
| void *ptr = NULL); | | | |
| void lsfitfit(lsfitstate &state, | | | |
| void (*func)(const real_1d_array &c, const real_1d_array &x, double &fu | | | |
| nc, void *ptr), | | | |
| void (*grad)(const real_1d_array &c, const real_1d_array &x, double &fu | | | |
| nc, real_1d_array &grad, void *ptr), | | | |
| void (*hess)(const real_1d_array &c, const real_1d_array &x, double &fu | | | |
| nc, real_1d_array &grad, real_2d_array &hess, void *ptr), | | | |
| void (*rep)(const real_1d_array &c, double func, void *ptr) = NULL, | | | |
| void *ptr = NULL); | | | |
| | | | |
| /************************************************************************* | | | |
| Nonlinear least squares fitting results. | | | |
| | | | |
| Called after return from LSFitFit(). | | | |
| | | | |
| INPUT PARAMETERS: | | | |
| State - algorithm state | | | |
| | | | |
| OUTPUT PARAMETERS: | | | |
| Info - completetion code: | | | |
| * 1 relative function improvement is no more than | | | |
| EpsF. | | | |
| * 2 relative step is no more than EpsX. | | | |
| * 4 gradient norm is no more than EpsG | | | |
| * 5 MaxIts steps was taken | | | |
| * 7 stopping conditions are too stringent, | | | |
| further improvement is impossible | | | |
| C - array[0..K-1], solution | | | |
| Rep - optimization report. Following fields are set: | | | |
| * Rep.TerminationType completetion code: | | | |
| * RMSError rms error on the (X,Y). | | | |
| * AvgError average error on the (X,Y). | | | |
| * AvgRelError average relative error on the non-zero | | | |
| Y | | | |
| * MaxError maximum error | | | |
| NON-WEIGHTED ERRORS ARE CALCULATED | | | |
| | | | |
| -- ALGLIB -- | | | |
| Copyright 17.08.2009 by Bochkanov Sergey | | | |
| *************************************************************************/ | | | |
| void lsfitresults(const lsfitstate &state, ae_int_t &info, real_1d_array &c | | | |
| , lsfitreport &rep); | | | |
| | | | |
| /************************************************************************* | | | |
| Lagrange intepolant: generation of the model on the general grid. | | | |
| This function has O(N^2) complexity. | | | |
| | | | |
| INPUT PARAMETERS: | | | |
| X - abscissas, array[0..N-1] | | | |
| Y - function values, array[0..N-1] | | | |
| N - number of points, N>=1 | | | |
| | | | |
| OIYTPUT PARAMETERS | | | |
| P - barycentric model which represents Lagrange interpolant | | | |
| (see ratint unit info and BarycentricCalc() description for | | | |
| more information). | | | |
| | | | |
| -- ALGLIB -- | | | |
| Copyright 02.12.2009 by Bochkanov Sergey | | | |
| *************************************************************************/ | | | |
| void polynomialbuild(const real_1d_array &x, const real_1d_array &y, const | | | |
| ae_int_t n, barycentricinterpolant &p); | | | |
| void polynomialbuild(const real_1d_array &x, const real_1d_array &y, baryce | | | |
| ntricinterpolant &p); | | | |
| | | | |
| /************************************************************************* | | | |
| Lagrange intepolant: generation of the model on equidistant grid. | | | |
| This function has O(N) complexity. | | | |
| | | | |
| INPUT PARAMETERS: | | | |
| A - left boundary of [A,B] | | | |
| B - right boundary of [A,B] | | | |
| Y - function values at the nodes, array[0..N-1] | | | |
| N - number of points, N>=1 | | | |
| for N=1 a constant model is constructed. | | | |
| | | | |
| OIYTPUT PARAMETERS | | | |
| P - barycentric model which represents Lagrange interpolant | | | |
| (see ratint unit info and BarycentricCalc() description for | | | |
| more information). | | | |
| | | | |
| -- ALGLIB -- | | | |
| Copyright 03.12.2009 by Bochkanov Sergey | | | |
| *************************************************************************/ | | | |
| void polynomialbuildeqdist(const double a, const double b, const real_1d_ar | | | |
| ray &y, const ae_int_t n, barycentricinterpolant &p); | | | |
| void polynomialbuildeqdist(const double a, const double b, const real_1d_ar | | | |
| ray &y, barycentricinterpolant &p); | | | |
| | | | |
| /************************************************************************* | | | |
| Lagrange intepolant on Chebyshev grid (first kind). | | | |
| This function has O(N) complexity. | | | |
| | | | |
| INPUT PARAMETERS: | | | |
| A - left boundary of [A,B] | | | |
| B - right boundary of [A,B] | | | |
| Y - function values at the nodes, array[0..N-1], | | | |
| Y[I] = Y(0.5*(B+A) + 0.5*(B-A)*Cos(PI*(2*i+1)/(2*n))) | | | |
| N - number of points, N>=1 | | | |
| for N=1 a constant model is constructed. | | | |
| | | | |
| OIYTPUT PARAMETERS | | | |
| P - barycentric model which represents Lagrange interpolant | | | |
| (see ratint unit info and BarycentricCalc() description for | | | |
| more information). | | | |
| | | | |
| -- ALGLIB -- | | | |
| Copyright 03.12.2009 by Bochkanov Sergey | | | |
| *************************************************************************/ | | | |
| void polynomialbuildcheb1(const double a, const double b, const real_1d_arr | | | |
| ay &y, const ae_int_t n, barycentricinterpolant &p); | | | |
| void polynomialbuildcheb1(const double a, const double b, const real_1d_arr | | | |
| ay &y, barycentricinterpolant &p); | | | |
| | | | |
| /************************************************************************* | | | |
| Lagrange intepolant on Chebyshev grid (second kind). | | | |
| This function has O(N) complexity. | | | |
| | | | |
| INPUT PARAMETERS: | | | |
| A - left boundary of [A,B] | | | |
| B - right boundary of [A,B] | | | |
| Y - function values at the nodes, array[0..N-1], | | | |
| Y[I] = Y(0.5*(B+A) + 0.5*(B-A)*Cos(PI*i/(n-1))) | | | |
| N - number of points, N>=1 | | | |
| for N=1 a constant model is constructed. | | | |
| | | | |
| OIYTPUT PARAMETERS | | | |
| P - barycentric model which represents Lagrange interpolant | | | |
| (see ratint unit info and BarycentricCalc() description for | | | |
| more information). | | | |
| | | | |
| -- ALGLIB -- | | | |
| Copyright 03.12.2009 by Bochkanov Sergey | | | |
| *************************************************************************/ | | | |
| void polynomialbuildcheb2(const double a, const double b, const real_1d_arr | | | |
| ay &y, const ae_int_t n, barycentricinterpolant &p); | | | |
| void polynomialbuildcheb2(const double a, const double b, const real_1d_arr | | | |
| ay &y, barycentricinterpolant &p); | | | |
| | | | |
| /************************************************************************* | | | |
| Fast equidistant polynomial interpolation function with O(N) complexity | | | |
| | | | |
| INPUT PARAMETERS: | | | |
| A - left boundary of [A,B] | | | |
| B - right boundary of [A,B] | | | |
| F - function values, array[0..N-1] | | | |
| N - number of points on equidistant grid, N>=1 | | | |
| for N=1 a constant model is constructed. | | | |
| T - position where P(x) is calculated | | | |
| | | | |
| RESULT | | | |
| value of the Lagrange interpolant at T | | | |
| | | | |
| IMPORTANT | | | |
| this function provides fast interface which is not overflow-safe | | | |
| nor it is very precise. | | | |
| the best option is to use PolynomialBuildEqDist()/BarycentricCalc() | | | |
| subroutines unless you are pretty sure that your data will not result | | | |
| in overflow. | | | |
| | | | |
| -- ALGLIB -- | | | |
| Copyright 02.12.2009 by Bochkanov Sergey | | | |
| *************************************************************************/ | | | |
| double polynomialcalceqdist(const double a, const double b, const real_1d_a | | | |
| rray &f, const ae_int_t n, const double t); | | | |
| double polynomialcalceqdist(const double a, const double b, const real_1d_a | | | |
| rray &f, const double t); | | | |
| | | | |
| /************************************************************************* | | | |
| Fast polynomial interpolation function on Chebyshev points (first kind) | | | |
| with O(N) complexity. | | | |
| | | | |
| INPUT PARAMETERS: | | | |
| A - left boundary of [A,B] | | | |
| B - right boundary of [A,B] | | | |
| F - function values, array[0..N-1] | | | |
| N - number of points on Chebyshev grid (first kind), | | | |
| X[i] = 0.5*(B+A) + 0.5*(B-A)*Cos(PI*(2*i+1)/(2*n)) | | | |
| for N=1 a constant model is constructed. | | | |
| T - position where P(x) is calculated | | | |
| | | | |
| RESULT | | | |
| value of the Lagrange interpolant at T | | | |
| | | | |
| IMPORTANT | | | |
| this function provides fast interface which is not overflow-safe | | | |
| nor it is very precise. | | | |
| the best option is to use PolIntBuildCheb1()/BarycentricCalc() | | | |
| subroutines unless you are pretty sure that your data will not result | | | |
| in overflow. | | | |
| | | | |
|
| -- ALGLIB -- | | For example, if we want to do linear fit by f(c0,c1,x) = c0*x+c1, then | |
| Copyright 02.12.2009 by Bochkanov Sergey | | x will be argument, and {c0,c1} will be parameters. | |
| *************************************************************************/ | | | |
| double polynomialcalccheb1(const double a, const double b, const real_1d_ar | | | |
| ray &f, const ae_int_t n, const double t); | | | |
| double polynomialcalccheb1(const double a, const double b, const real_1d_ar | | | |
| ray &f, const double t); | | | |
| | | | |
|
| /************************************************************************* | | It is important to understand that this algorithm finds minimum in the | |
| Fast polynomial interpolation function on Chebyshev points (second kind) | | space of function PARAMETERS (not arguments), so it needs derivatives | |
| with O(N) complexity. | | of f() with respect to C, not X. | |
| | | | |
|
| INPUT PARAMETERS: | | In the example above it will need f=c0*x+c1 and {df/dc0,df/dc1} = {x,1} | |
| A - left boundary of [A,B] | | instead of {df/dx} = {c0}. | |
| B - right boundary of [A,B] | | | |
| F - function values, array[0..N-1] | | | |
| N - number of points on Chebyshev grid (second kind), | | | |
| X[i] = 0.5*(B+A) + 0.5*(B-A)*Cos(PI*i/(n-1)) | | | |
| for N=1 a constant model is constructed. | | | |
| T - position where P(x) is calculated | | | |
| | | | |
|
| RESULT | | 2. Callback functions accept C as the first parameter, and X as the second | |
| value of the Lagrange interpolant at T | | | |
| | | | |
|
| IMPORTANT | | 3. If state was created with LSFitCreateFG(), algorithm needs just | |
| this function provides fast interface which is not overflow-safe | | function and its gradient, but if state was created with | |
| nor it is very precise. | | LSFitCreateFGH(), algorithm will need function, gradient and Hessian. | |
| the best option is to use PolIntBuildCheb2()/BarycentricCalc() | | | |
| subroutines unless you are pretty sure that your data will not result | | | |
| in overflow. | | | |
| | | | |
|
| -- ALGLIB -- | | According to the said above, there ase several versions of this | |
| Copyright 02.12.2009 by Bochkanov Sergey | | function, which accept different sets of callbacks. | |
| *************************************************************************/ | | | |
| double polynomialcalccheb2(const double a, const double b, const real_1d_ar | | | |
| ray &f, const ae_int_t n, const double t); | | | |
| double polynomialcalccheb2(const double a, const double b, const real_1d_ar | | | |
| ray &f, const double t); | | | |
| | | | |
|
| /************************************************************************* | | This flexibility opens way to subtle errors - you may create state with | |
| Least squares fitting by polynomial. | | LSFitCreateFGH() (optimization using Hessian), but call function which | |
| | | does not accept Hessian. So when algorithm will request Hessian, there | |
| | | will be no callback to call. In this case exception will be thrown. | |
| | | | |
|
| This subroutine is "lightweight" alternative for more complex and feature- | | Be careful to avoid such errors because there is no way to find them at | |
| rich PolynomialFitWC(). See PolynomialFitWC() for more information about | | compile time - you can see them at runtime only. | |
| subroutine parameters (we don't duplicate it here because of length) | | | |
| | | -- ALGLIB -- | |
| | | Copyright 17.08.2009 by Bochkanov Sergey | |
| | | | |
|
| -- ALGLIB PROJECT -- | | | |
| Copyright 12.10.2009 by Bochkanov Sergey | | | |
| *************************************************************************/ | | *************************************************************************/ | |
|
| void polynomialfit(const real_1d_array &x, const real_1d_array &y, const ae | | void lsfitfit(lsfitstate &state, | |
| _int_t n, const ae_int_t m, ae_int_t &info, barycentricinterpolant &p, poly | | void (*func)(const real_1d_array &c, const real_1d_array &x, double &fu | |
| nomialfitreport &rep); | | nc, void *ptr), | |
| void polynomialfit(const real_1d_array &x, const real_1d_array &y, const ae | | void (*rep)(const real_1d_array &c, double func, void *ptr) = NULL, | |
| _int_t m, ae_int_t &info, barycentricinterpolant &p, polynomialfitreport &r | | void *ptr = NULL); | |
| ep); | | void lsfitfit(lsfitstate &state, | |
| | | void (*func)(const real_1d_array &c, const real_1d_array &x, double &fu | |
| | | nc, void *ptr), | |
| | | void (*grad)(const real_1d_array &c, const real_1d_array &x, double &fu | |
| | | nc, real_1d_array &grad, void *ptr), | |
| | | void (*rep)(const real_1d_array &c, double func, void *ptr) = NULL, | |
| | | void *ptr = NULL); | |
| | | void lsfitfit(lsfitstate &state, | |
| | | void (*func)(const real_1d_array &c, const real_1d_array &x, double &fu | |
| | | nc, void *ptr), | |
| | | void (*grad)(const real_1d_array &c, const real_1d_array &x, double &fu | |
| | | nc, real_1d_array &grad, void *ptr), | |
| | | void (*hess)(const real_1d_array &c, const real_1d_array &x, double &fu | |
| | | nc, real_1d_array &grad, real_2d_array &hess, void *ptr), | |
| | | void (*rep)(const real_1d_array &c, double func, void *ptr) = NULL, | |
| | | void *ptr = NULL); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
|
| Weighted fitting by Chebyshev polynomial in barycentric form, with | | Nonlinear least squares fitting results. | |
| constraints on function values or first derivatives. | | | |
| | | | |
| Small regularizing term is used when solving constrained tasks (to improve | | | |
| stability). | | | |
| | | | |
| Task is linear, so linear least squares solver is used. Complexity of this | | | |
| computational scheme is O(N*M^2), mostly dominated by least squares solver | | | |
| | | | |
|
| SEE ALSO: | | Called after return from LSFitFit(). | |
| PolynomialFit() | | | |
| | | | |
| INPUT PARAMETERS: | | INPUT PARAMETERS: | |
|
| X - points, array[0..N-1]. | | State - algorithm state | |
| Y - function values, array[0..N-1]. | | | |
| W - weights, array[0..N-1] | | | |
| Each summand in square sum of approximation deviations from | | | |
| given values is multiplied by the square of corresponding | | | |
| weight. Fill it by 1's if you don't want to solve weighted | | | |
| task. | | | |
| N - number of points, N>0. | | | |
| XC - points where polynomial values/derivatives are constrained, | | | |
| array[0..K-1]. | | | |
| YC - values of constraints, array[0..K-1] | | | |
| DC - array[0..K-1], types of constraints: | | | |
| * DC[i]=0 means that P(XC[i])=YC[i] | | | |
| * DC[i]=1 means that P'(XC[i])=YC[i] | | | |
| SEE BELOW FOR IMPORTANT INFORMATION ON CONSTRAINTS | | | |
| K - number of constraints, 0<=K<M. | | | |
| K=0 means no constraints (XC/YC/DC are not used in such cases) | | | |
| M - number of basis functions (= polynomial_degree + 1), M>=1 | | | |
| | | | |
| OUTPUT PARAMETERS: | | OUTPUT PARAMETERS: | |
|
| Info- same format as in LSFitLinearW() subroutine: | | Info - completetion code: | |
| * Info>0 task is solved | | * 1 relative function improvement is no more than | |
| * Info<=0 an error occured: | | EpsF. | |
| -4 means inconvergence of internal SVD | | * 2 relative step is no more than EpsX. | |
| -3 means inconsistent constraints | | * 4 gradient norm is no more than EpsG | |
| P - interpolant in barycentric form. | | * 5 MaxIts steps was taken | |
| Rep - report, same format as in LSFitLinearW() subroutine. | | * 7 stopping conditions are too stringent, | |
| Following fields are set: | | further improvement is impossible | |
| * RMSError rms error on the (X,Y). | | C - array[0..K-1], solution | |
| * AvgError average error on the (X,Y). | | Rep - optimization report. Following fields are set: | |
| * AvgRelError average relative error on the non-zero Y | | * Rep.TerminationType completetion code: | |
| * MaxError maximum error | | * RMSError rms error on the (X,Y). | |
| NON-WEIGHTED ERRORS ARE CALCULATED | | * AvgError average error on the (X,Y). | |
| | | * AvgRelError average relative error on the non-zero | |
| IMPORTANT: | | Y | |
| this subroitine doesn't calculate task's condition number for K<>0. | | * MaxError maximum error | |
| | | NON-WEIGHTED ERRORS ARE CALCULATED | |
| SETTING CONSTRAINTS - DANGERS AND OPPORTUNITIES: | | | |
| | | | |
| Setting constraints can lead to undesired results, like ill-conditioned | | | |
| behavior, or inconsistency being detected. From the other side, it allows | | | |
| us to improve quality of the fit. Here we summarize our experience with | | | |
| constrained regression splines: | | | |
| * even simple constraints can be inconsistent, see Wikipedia article on | | | |
| this subject: http://en.wikipedia.org/wiki/Birkhoff_interpolation | | | |
| * the greater is M (given fixed constraints), the more chances that | | | |
| constraints will be consistent | | | |
| * in the general case, consistency of constraints is NOT GUARANTEED. | | | |
| * in the one special cases, however, we can guarantee consistency. This | | | |
| case is: M>1 and constraints on the function values (NOT DERIVATIVES) | | | |
| | | | |
| Our final recommendation is to use constraints WHEN AND ONLY when you | | | |
| can't solve your task without them. Anything beyond special cases given | | | |
| above is not guaranteed and may result in inconsistency. | | | |
| | | | |
|
| -- ALGLIB PROJECT -- | | -- ALGLIB -- | |
| Copyright 10.12.2009 by Bochkanov Sergey | | Copyright 17.08.2009 by Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
|
| void polynomialfitwc(const real_1d_array &x, const real_1d_array &y, const | | void lsfitresults(const lsfitstate &state, ae_int_t &info, real_1d_array &c | |
| real_1d_array &w, const ae_int_t n, const real_1d_array &xc, const real_1d_ | | , lsfitreport &rep); | |
| array &yc, const integer_1d_array &dc, const ae_int_t k, const ae_int_t m, | | | |
| ae_int_t &info, barycentricinterpolant &p, polynomialfitreport &rep); | | | |
| void polynomialfitwc(const real_1d_array &x, const real_1d_array &y, const | | | |
| real_1d_array &w, const real_1d_array &xc, const real_1d_array &yc, const i | | | |
| nteger_1d_array &dc, const ae_int_t m, ae_int_t &info, barycentricinterpola | | | |
| nt &p, polynomialfitreport &rep); | | | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
| Rational interpolation using barycentric formula | | Rational interpolation using barycentric formula | |
| | | | |
| F(t) = SUM(i=0,n-1,w[i]*f[i]/(t-x[i])) / SUM(i=0,n-1,w[i]/(t-x[i])) | | F(t) = SUM(i=0,n-1,w[i]*f[i]/(t-x[i])) / SUM(i=0,n-1,w[i]/(t-x[i])) | |
| | | | |
| Input parameters: | | Input parameters: | |
| B - barycentric interpolant built with one of model building | | B - barycentric interpolant built with one of model building | |
| subroutines. | | subroutines. | |
| T - interpolation point | | T - interpolation point | |
| | | | |
| skipping to change at line 1364 | | skipping to change at line 1829 | |
|---|
| Note: | | Note: | |
| this algorithm always succeeds and calculates the weights with close | | this algorithm always succeeds and calculates the weights with close | |
| to machine precision. | | to machine precision. | |
| | | | |
| -- ALGLIB PROJECT -- | | -- ALGLIB PROJECT -- | |
| Copyright 17.06.2007 by Bochkanov Sergey | | Copyright 17.06.2007 by Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
| void barycentricbuildfloaterhormann(const real_1d_array &x, const real_1d_a
rray &y, const ae_int_t n, const ae_int_t d, barycentricinterpolant &b); | | void barycentricbuildfloaterhormann(const real_1d_array &x, const real_1d_a
rray &y, const ae_int_t n, const ae_int_t d, barycentricinterpolant &b); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
|
| Weghted rational least squares fitting using Floater-Hormann rational | | Conversion from barycentric representation to Chebyshev basis. | |
| functions with optimal D chosen from [0,9], with constraints and | | This function has O(N^2) complexity. | |
| individual weights. | | | |
| | | | |
| Equidistant grid with M node on [min(x),max(x)] is used to build basis | | | |
| functions. Different values of D are tried, optimal D (least WEIGHTED root | | | |
| mean square error) is chosen. Task is linear, so linear least squares | | | |
| solver is used. Complexity of this computational scheme is O(N*M^2) | | | |
| (mostly dominated by the least squares solver). | | | |
| | | | |
| SEE ALSO | | | |
| * BarycentricFitFloaterHormann(), "lightweight" fitting without invididual | | | |
| weights and constraints. | | | |
| | | | |
| INPUT PARAMETERS: | | INPUT PARAMETERS: | |
|
| X - points, array[0..N-1]. | | P - polynomial in barycentric form | |
| Y - function values, array[0..N-1]. | | A,B - base interval for Chebyshev polynomials (see below) | |
| W - weights, array[0..N-1] | | A<>B | |
| Each summand in square sum of approximation deviations from | | | |
| given values is multiplied by the square of corresponding | | OUTPUT PARAMETERS | |
| weight. Fill it by 1's if you don't want to solve weighted | | T - coefficients of Chebyshev representation; | |
| task. | | P(x) = sum { T[i]*Ti(2*(x-A)/(B-A)-1), i=0..N-1 }, | |
| N - number of points, N>0. | | where Ti - I-th Chebyshev polynomial. | |
| XC - points where function values/derivatives are constrained, | | | |
| array[0..K-1]. | | | |
| YC - values of constraints, array[0..K-1] | | | |
| DC - array[0..K-1], types of constraints: | | | |
| * DC[i]=0 means that S(XC[i])=YC[i] | | | |
| * DC[i]=1 means that S'(XC[i])=YC[i] | | | |
| SEE BELOW FOR IMPORTANT INFORMATION ON CONSTRAINTS | | | |
| K - number of constraints, 0<=K<M. | | | |
| K=0 means no constraints (XC/YC/DC are not used in such cases) | | | |
| M - number of basis functions ( = number_of_nodes), M>=2. | | | |
| | | | |
|
| OUTPUT PARAMETERS: | | NOTES: | |
| Info- same format as in LSFitLinearWC() subroutine. | | barycentric interpolant passed as P may be either polynomial obtained | |
| * Info>0 task is solved | | from polynomial interpolation/ fitting or rational function which is | |
| * Info<=0 an error occured: | | NOT polynomial. We can't distinguish between these two cases, and this | |
| -4 means inconvergence of internal SVD | | algorithm just tries to work assuming that P IS a polynomial. If not, | |
| -3 means inconsistent constraints | | algorithm will return results, but they won't have any meaning. | |
| -1 means another errors in parameters passed | | | |
| (N<=0, for example) | | | |
| B - barycentric interpolant. | | | |
| Rep - report, same format as in LSFitLinearWC() subroutine. | | | |
| Following fields are set: | | | |
| * DBest best value of the D parameter | | | |
| * RMSError rms error on the (X,Y). | | | |
| * AvgError average error on the (X,Y). | | | |
| * AvgRelError average relative error on the non-zero Y | | | |
| * MaxError maximum error | | | |
| NON-WEIGHTED ERRORS ARE CALCULATED | | | |
| | | | |
|
| IMPORTANT: | | -- ALGLIB -- | |
| this subroitine doesn't calculate task's condition number for K<>0. | | Copyright 30.09.2010 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void polynomialbar2cheb(const barycentricinterpolant &p, const double a, co | |
| | | nst double b, real_1d_array &t); | |
| | | | |
|
| SETTING CONSTRAINTS - DANGERS AND OPPORTUNITIES: | | /************************************************************************* | |
| | | Conversion from Chebyshev basis to barycentric representation. | |
| | | This function has O(N^2) complexity. | |
| | | | |
|
| Setting constraints can lead to undesired results, like ill-conditioned | | INPUT PARAMETERS: | |
| behavior, or inconsistency being detected. From the other side, it allows | | T - coefficients of Chebyshev representation; | |
| us to improve quality of the fit. Here we summarize our experience with | | P(x) = sum { T[i]*Ti(2*(x-A)/(B-A)-1), i=0..N }, | |
| constrained barycentric interpolants: | | where Ti - I-th Chebyshev polynomial. | |
| * excessive constraints can be inconsistent. Floater-Hormann basis | | N - number of coefficients: | |
| functions aren't as flexible as splines (although they are very smooth). | | * if given, only leading N elements of T are used | |
| * the more evenly constraints are spread across [min(x),max(x)], the more | | * if not given, automatically determined from size of T | |
| chances that they will be consistent | | A,B - base interval for Chebyshev polynomials (see above) | |
| * the greater is M (given fixed constraints), the more chances that | | A<B | |
| constraints will be consistent | | | |
| * in the general case, consistency of constraints IS NOT GUARANTEED. | | | |
| * in the several special cases, however, we CAN guarantee consistency. | | | |
| * one of this cases is constraints on the function VALUES at the interval | | | |
| boundaries. Note that consustency of the constraints on the function | | | |
| DERIVATIVES is NOT guaranteed (you can use in such cases cubic splines | | | |
| which are more flexible). | | | |
| * another special case is ONE constraint on the function value (OR, but | | | |
| not AND, derivative) anywhere in the interval | | | |
| | | | |
|
| Our final recommendation is to use constraints WHEN AND ONLY WHEN you | | OUTPUT PARAMETERS | |
| can't solve your task without them. Anything beyond special cases given | | P - polynomial in barycentric form | |
| above is not guaranteed and may result in inconsistency. | | | |
| | | | |
|
| -- ALGLIB PROJECT -- | | -- ALGLIB -- | |
| Copyright 18.08.2009 by Bochkanov Sergey | | Copyright 30.09.2010 by Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
|
| void barycentricfitfloaterhormannwc(const real_1d_array &x, const real_1d_a | | void polynomialcheb2bar(const real_1d_array &t, const ae_int_t n, const dou | |
| rray &y, const real_1d_array &w, const ae_int_t n, const real_1d_array &xc, | | ble a, const double b, barycentricinterpolant &p); | |
| const real_1d_array &yc, const integer_1d_array &dc, const ae_int_t k, con | | void polynomialcheb2bar(const real_1d_array &t, const double a, const doubl | |
| st ae_int_t m, ae_int_t &info, barycentricinterpolant &b, barycentricfitrep | | e b, barycentricinterpolant &p); | |
| ort &rep); | | | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
|
| Rational least squares fitting, without weights and constraints. | | Conversion from barycentric representation to power basis. | |
| | | This function has O(N^2) complexity. | |
| | | | |
|
| See BarycentricFitFloaterHormannWC() for more information. | | INPUT PARAMETERS: | |
| | | P - polynomial in barycentric form | |
| | | C - offset (see below); 0.0 is used as default value. | |
| | | S - scale (see below); 1.0 is used as default value. S<>0. | |
| | | | |
| | | OUTPUT PARAMETERS | |
| | | A - coefficients, P(x) = sum { A[i]*((X-C)/S)^i, i=0..N-1 } | |
| | | N - number of coefficients (polynomial degree plus 1) | |
| | | | |
|
| -- ALGLIB PROJECT -- | | NOTES: | |
| Copyright 18.08.2009 by Bochkanov Sergey | | 1. this function accepts offset and scale, which can be set to improve | |
| *************************************************************************/ | | numerical properties of polynomial. For example, if P was obtained as | |
| void barycentricfitfloaterhormann(const real_1d_array &x, const real_1d_arr | | result of interpolation on [-1,+1], you can set C=0 and S=1 and | |
| ay &y, const ae_int_t n, const ae_int_t m, ae_int_t &info, barycentricinter | | represent P as sum of 1, x, x^2, x^3 and so on. In most cases you it | |
| polant &b, barycentricfitreport &rep); | | is exactly what you need. | |
| | | | |
|
| /************************************************************************* | | However, if your interpolation model was built on [999,1001], you will | |
| This subroutine builds bilinear spline coefficients table. | | see significant growth of numerical errors when using {1, x, x^2, x^3} | |
| | | as basis. Representing P as sum of 1, (x-1000), (x-1000)^2, (x-1000)^3 | |
| | | will be better option. Such representation can be obtained by using | |
| | | 1000.0 as offset C and 1.0 as scale S. | |
| | | | |
|
| Input parameters: | | 2. power basis is ill-conditioned and tricks described above can't solve | |
| X - spline abscissas, array[0..N-1] | | this problem completely. This function will return coefficients in | |
| Y - spline ordinates, array[0..M-1] | | any case, but for N>8 they will become unreliable. However, N's | |
| F - function values, array[0..M-1,0..N-1] | | less than 5 are pretty safe. | |
| M,N - grid size, M>=2, N>=2 | | | |
| | | | |
|
| Output parameters: | | 3. barycentric interpolant passed as P may be either polynomial obtained | |
| C - spline interpolant | | from polynomial interpolation/ fitting or rational function which is | |
| | | NOT polynomial. We can't distinguish between these two cases, and this | |
| | | algorithm just tries to work assuming that P IS a polynomial. If not, | |
| | | algorithm will return results, but they won't have any meaning. | |
| | | | |
|
| -- ALGLIB PROJECT -- | | -- ALGLIB -- | |
| Copyright 05.07.2007 by Bochkanov Sergey | | Copyright 30.09.2010 by Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
|
| void spline2dbuildbilinear(const real_1d_array &x, const real_1d_array &y, | | void polynomialbar2pow(const barycentricinterpolant &p, const double c, con | |
| const real_2d_array &f, const ae_int_t m, const ae_int_t n, spline2dinterpo | | st double s, real_1d_array &a); | |
| lant &c); | | void polynomialbar2pow(const barycentricinterpolant &p, real_1d_array &a); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
|
| This subroutine builds bicubic spline coefficients table. | | Conversion from power basis to barycentric representation. | |
| | | This function has O(N^2) complexity. | |
| | | | |
|
| Input parameters: | | INPUT PARAMETERS: | |
| X - spline abscissas, array[0..N-1] | | A - coefficients, P(x) = sum { A[i]*((X-C)/S)^i, i=0..N-1 } | |
| Y - spline ordinates, array[0..M-1] | | N - number of coefficients (polynomial degree plus 1) | |
| F - function values, array[0..M-1,0..N-1] | | * if given, only leading N elements of A are used | |
| M,N - grid size, M>=2, N>=2 | | * if not given, automatically determined from size of A | |
| | | C - offset (see below); 0.0 is used as default value. | |
| | | S - scale (see below); 1.0 is used as default value. S<>0. | |
| | | | |
|
| Output parameters: | | OUTPUT PARAMETERS | |
| C - spline interpolant | | P - polynomial in barycentric form | |
| | | | |
|
| -- ALGLIB PROJECT -- | | NOTES: | |
| Copyright 05.07.2007 by Bochkanov Sergey | | 1. this function accepts offset and scale, which can be set to improve | |
| | | numerical properties of polynomial. For example, if you interpolate on | |
| | | [-1,+1], you can set C=0 and S=1 and convert from sum of 1, x, x^2, | |
| | | x^3 and so on. In most cases you it is exactly what you need. | |
| | | | |
| | | However, if your interpolation model was built on [999,1001], you will | |
| | | see significant growth of numerical errors when using {1, x, x^2, x^3} | |
| | | as input basis. Converting from sum of 1, (x-1000), (x-1000)^2, | |
| | | (x-1000)^3 will be better option (you have to specify 1000.0 as offset | |
| | | C and 1.0 as scale S). | |
| | | | |
| | | 2. power basis is ill-conditioned and tricks described above can't solve | |
| | | this problem completely. This function will return barycentric model | |
| | | in any case, but for N>8 accuracy well degrade. However, N's less than | |
| | | 5 are pretty safe. | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 30.09.2010 by Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
|
| void spline2dbuildbicubic(const real_1d_array &x, const real_1d_array &y, c | | void polynomialpow2bar(const real_1d_array &a, const ae_int_t n, const doub | |
| onst real_2d_array &f, const ae_int_t m, const ae_int_t n, spline2dinterpol | | le c, const double s, barycentricinterpolant &p); | |
| ant &c); | | void polynomialpow2bar(const real_1d_array &a, barycentricinterpolant &p); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
|
| This subroutine calculates the value of the bilinear or bicubic spline at | | Lagrange intepolant: generation of the model on the general grid. | |
| the given point X. | | This function has O(N^2) complexity. | |
| | | | |
|
| Input parameters: | | INPUT PARAMETERS: | |
| C - coefficients table. | | X - abscissas, array[0..N-1] | |
| Built by BuildBilinearSpline or BuildBicubicSpline. | | Y - function values, array[0..N-1] | |
| X, Y- point | | N - number of points, N>=1 | |
| | | | |
|
| Result: | | OUTPUT PARAMETERS | |
| S(x,y) | | P - barycentric model which represents Lagrange interpolant | |
| | | (see ratint unit info and BarycentricCalc() description for | |
| | | more information). | |
| | | | |
|
| -- ALGLIB PROJECT -- | | -- ALGLIB -- | |
| Copyright 05.07.2007 by Bochkanov Sergey | | Copyright 02.12.2009 by Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
|
| double spline2dcalc(const spline2dinterpolant &c, const double x, const dou | | void polynomialbuild(const real_1d_array &x, const real_1d_array &y, const | |
| ble y); | | ae_int_t n, barycentricinterpolant &p); | |
| | | void polynomialbuild(const real_1d_array &x, const real_1d_array &y, baryce | |
| | | ntricinterpolant &p); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
|
| This subroutine calculates the value of the bilinear or bicubic spline at | | Lagrange intepolant: generation of the model on equidistant grid. | |
| the given point X and its derivatives. | | This function has O(N) complexity. | |
| | | | |
|
| Input parameters: | | INPUT PARAMETERS: | |
| C - spline interpolant. | | A - left boundary of [A,B] | |
| X, Y- point | | B - right boundary of [A,B] | |
| | | Y - function values at the nodes, array[0..N-1] | |
| | | N - number of points, N>=1 | |
| | | for N=1 a constant model is constructed. | |
| | | | |
|
| Output parameters: | | OUTPUT PARAMETERS | |
| F - S(x,y) | | P - barycentric model which represents Lagrange interpolant | |
| FX - dS(x,y)/dX | | (see ratint unit info and BarycentricCalc() description for | |
| FY - dS(x,y)/dY | | more information). | |
| FXY - d2S(x,y)/dXdY | | | |
| | | | |
|
| -- ALGLIB PROJECT -- | | -- ALGLIB -- | |
| Copyright 05.07.2007 by Bochkanov Sergey | | Copyright 03.12.2009 by Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
|
| void spline2ddiff(const spline2dinterpolant &c, const double x, const doubl | | void polynomialbuildeqdist(const double a, const double b, const real_1d_ar | |
| e y, double &f, double &fx, double &fy, double &fxy); | | ray &y, const ae_int_t n, barycentricinterpolant &p); | |
| | | void polynomialbuildeqdist(const double a, const double b, const real_1d_ar | |
| | | ray &y, barycentricinterpolant &p); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
|
| This subroutine unpacks two-dimensional spline into the coefficients table | | Lagrange intepolant on Chebyshev grid (first kind). | |
| | | This function has O(N) complexity. | |
| | | | |
|
| Input parameters: | | INPUT PARAMETERS: | |
| C - spline interpolant. | | A - left boundary of [A,B] | |
| | | B - right boundary of [A,B] | |
| | | Y - function values at the nodes, array[0..N-1], | |
| | | Y[I] = Y(0.5*(B+A) + 0.5*(B-A)*Cos(PI*(2*i+1)/(2*n))) | |
| | | N - number of points, N>=1 | |
| | | for N=1 a constant model is constructed. | |
| | | | |
|
| Result: | | OUTPUT PARAMETERS | |
| M, N- grid size (x-axis and y-axis) | | P - barycentric model which represents Lagrange interpolant | |
| Tbl - coefficients table, unpacked format, | | (see ratint unit info and BarycentricCalc() description for | |
| [0..(N-1)*(M-1)-1, 0..19]. | | more information). | |
| For I = 0...M-2, J=0..N-2: | | | |
| K = I*(N-1)+J | | | |
| Tbl[K,0] = X[j] | | | |
| Tbl[K,1] = X[j+1] | | | |
| Tbl[K,2] = Y[i] | | | |
| Tbl[K,3] = Y[i+1] | | | |
| Tbl[K,4] = C00 | | | |
| Tbl[K,5] = C01 | | | |
| Tbl[K,6] = C02 | | | |
| Tbl[K,7] = C03 | | | |
| Tbl[K,8] = C10 | | | |
| Tbl[K,9] = C11 | | | |
| ... | | | |
| Tbl[K,19] = C33 | | | |
| On each grid square spline is equals to: | | | |
| S(x) = SUM(c[i,j]*(x^i)*(y^j), i=0..3, j=0..3) | | | |
| t = x-x[j] | | | |
| u = y-y[i] | | | |
| | | | |
|
| -- ALGLIB PROJECT -- | | -- ALGLIB -- | |
| Copyright 29.06.2007 by Bochkanov Sergey | | Copyright 03.12.2009 by Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
|
| void spline2dunpack(const spline2dinterpolant &c, ae_int_t &m, ae_int_t &n, | | void polynomialbuildcheb1(const double a, const double b, const real_1d_arr | |
| real_2d_array &tbl); | | ay &y, const ae_int_t n, barycentricinterpolant &p); | |
| | | void polynomialbuildcheb1(const double a, const double b, const real_1d_arr | |
| | | ay &y, barycentricinterpolant &p); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
|
| This subroutine performs linear transformation of the spline argument. | | Lagrange intepolant on Chebyshev grid (second kind). | |
| | | This function has O(N) complexity. | |
| | | | |
|
| Input parameters: | | INPUT PARAMETERS: | |
| C - spline interpolant | | A - left boundary of [A,B] | |
| AX, BX - transformation coefficients: x = A*t + B | | B - right boundary of [A,B] | |
| AY, BY - transformation coefficients: y = A*u + B | | Y - function values at the nodes, array[0..N-1], | |
| Result: | | Y[I] = Y(0.5*(B+A) + 0.5*(B-A)*Cos(PI*i/(n-1))) | |
| C - transformed spline | | N - number of points, N>=1 | |
| | | for N=1 a constant model is constructed. | |
| | | | |
|
| -- ALGLIB PROJECT -- | | OUTPUT PARAMETERS | |
| Copyright 30.06.2007 by Bochkanov Sergey | | P - barycentric model which represents Lagrange interpolant | |
| | | (see ratint unit info and BarycentricCalc() description for | |
| | | more information). | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 03.12.2009 by Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
|
| void spline2dlintransxy(const spline2dinterpolant &c, const double ax, cons | | void polynomialbuildcheb2(const double a, const double b, const real_1d_arr | |
| t double bx, const double ay, const double by); | | ay &y, const ae_int_t n, barycentricinterpolant &p); | |
| | | void polynomialbuildcheb2(const double a, const double b, const real_1d_arr | |
| | | ay &y, barycentricinterpolant &p); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
|
| This subroutine performs linear transformation of the spline. | | Fast equidistant polynomial interpolation function with O(N) complexity | |
| | | | |
|
| Input parameters: | | INPUT PARAMETERS: | |
| C - spline interpolant. | | A - left boundary of [A,B] | |
| A, B- transformation coefficients: S2(x,y) = A*S(x,y) + B | | B - right boundary of [A,B] | |
| | | F - function values, array[0..N-1] | |
| | | N - number of points on equidistant grid, N>=1 | |
| | | for N=1 a constant model is constructed. | |
| | | T - position where P(x) is calculated | |
| | | | |
|
| Output parameters: | | RESULT | |
| C - transformed spline | | value of the Lagrange interpolant at T | |
| | | | |
|
| -- ALGLIB PROJECT -- | | IMPORTANT | |
| Copyright 30.06.2007 by Bochkanov Sergey | | this function provides fast interface which is not overflow-safe | |
| | | nor it is very precise. | |
| | | the best option is to use PolynomialBuildEqDist()/BarycentricCalc() | |
| | | subroutines unless you are pretty sure that your data will not result | |
| | | in overflow. | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 02.12.2009 by Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
|
| void spline2dlintransf(const spline2dinterpolant &c, const double a, const | | double polynomialcalceqdist(const double a, const double b, const real_1d_a | |
| double b); | | rray &f, const ae_int_t n, const double t); | |
| | | double polynomialcalceqdist(const double a, const double b, const real_1d_a | |
| | | rray &f, const double t); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
|
| Bicubic spline resampling | | Fast polynomial interpolation function on Chebyshev points (first kind) | |
| | | with O(N) complexity. | |
| | | | |
|
| Input parameters: | | INPUT PARAMETERS: | |
| A - function values at the old grid, | | A - left boundary of [A,B] | |
| array[0..OldHeight-1, 0..OldWidth-1] | | B - right boundary of [A,B] | |
| OldHeight - old grid height, OldHeight>1 | | F - function values, array[0..N-1] | |
| OldWidth - old grid width, OldWidth>1 | | N - number of points on Chebyshev grid (first kind), | |
| NewHeight - new grid height, NewHeight>1 | | X[i] = 0.5*(B+A) + 0.5*(B-A)*Cos(PI*(2*i+1)/(2*n)) | |
| NewWidth - new grid width, NewWidth>1 | | for N=1 a constant model is constructed. | |
| | | T - position where P(x) is calculated | |
| | | | |
|
| Output parameters: | | RESULT | |
| B - function values at the new grid, | | value of the Lagrange interpolant at T | |
| array[0..NewHeight-1, 0..NewWidth-1] | | | |
| | | | |
|
| -- ALGLIB routine -- | | IMPORTANT | |
| 15 May, 2007 | | this function provides fast interface which is not overflow-safe | |
| Copyright by Bochkanov Sergey | | nor it is very precise. | |
| | | the best option is to use PolIntBuildCheb1()/BarycentricCalc() | |
| | | subroutines unless you are pretty sure that your data will not result | |
| | | in overflow. | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 02.12.2009 by Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
|
| void spline2dresamplebicubic(const real_2d_array &a, const ae_int_t oldheig | | double polynomialcalccheb1(const double a, const double b, const real_1d_ar | |
| ht, const ae_int_t oldwidth, real_2d_array &b, const ae_int_t newheight, co | | ray &f, const ae_int_t n, const double t); | |
| nst ae_int_t newwidth); | | double polynomialcalccheb1(const double a, const double b, const real_1d_ar | |
| | | ray &f, const double t); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
|
| Bilinear spline resampling | | Fast polynomial interpolation function on Chebyshev points (second kind) | |
| | | with O(N) complexity. | |
| | | | |
|
| Input parameters: | | INPUT PARAMETERS: | |
| A - function values at the old grid, | | A - left boundary of [A,B] | |
| array[0..OldHeight-1, 0..OldWidth-1] | | B - right boundary of [A,B] | |
| OldHeight - old grid height, OldHeight>1 | | F - function values, array[0..N-1] | |
| OldWidth - old grid width, OldWidth>1 | | N - number of points on Chebyshev grid (second kind), | |
| NewHeight - new grid height, NewHeight>1 | | X[i] = 0.5*(B+A) + 0.5*(B-A)*Cos(PI*i/(n-1)) | |
| NewWidth - new grid width, NewWidth>1 | | for N=1 a constant model is constructed. | |
| | | T - position where P(x) is calculated | |
| | | | |
|
| Output parameters: | | RESULT | |
| B - function values at the new grid, | | value of the Lagrange interpolant at T | |
| array[0..NewHeight-1, 0..NewWidth-1] | | | |
| | | | |
|
| -- ALGLIB routine -- | | IMPORTANT | |
| 09.07.2007 | | this function provides fast interface which is not overflow-safe | |
| Copyright by Bochkanov Sergey | | nor it is very precise. | |
| | | the best option is to use PolIntBuildCheb2()/BarycentricCalc() | |
| | | subroutines unless you are pretty sure that your data will not result | |
| | | in overflow. | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 02.12.2009 by Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
|
| void spline2dresamplebilinear(const real_2d_array &a, const ae_int_t oldhei | | double polynomialcalccheb2(const double a, const double b, const real_1d_ar | |
| ght, const ae_int_t oldwidth, real_2d_array &b, const ae_int_t newheight, c | | ray &f, const ae_int_t n, const double t); | |
| onst ae_int_t newwidth); | | double polynomialcalccheb2(const double a, const double b, const real_1d_ar | |
| | | ray &f, const double t); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
| This subroutine builds linear spline interpolant | | This subroutine builds linear spline interpolant | |
| | | | |
| INPUT PARAMETERS: | | INPUT PARAMETERS: | |
| X - spline nodes, array[0..N-1] | | X - spline nodes, array[0..N-1] | |
| Y - function values, array[0..N-1] | | Y - function values, array[0..N-1] | |
| N - points count (optional): | | N - points count (optional): | |
| * N>=2 | | * N>=2 | |
| * if given, only first N points are used to build spline | | * if given, only first N points are used to build spline | |
| | | | |
| skipping to change at line 2059 | | skipping to change at line 2555 | |
|---|
| | | | |
| -- ALGLIB PROJECT -- | | -- ALGLIB PROJECT -- | |
| Copyright 03.09.2010 by Bochkanov Sergey | | Copyright 03.09.2010 by Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
| void spline1dconvdiff2cubic(const real_1d_array &x, const real_1d_array &y,
const ae_int_t n, const ae_int_t boundltype, const double boundl, const ae
_int_t boundrtype, const double boundr, const real_1d_array &x2, const ae_i
nt_t n2, real_1d_array &y2, real_1d_array &d2, real_1d_array &dd2); | | void spline1dconvdiff2cubic(const real_1d_array &x, const real_1d_array &y,
const ae_int_t n, const ae_int_t boundltype, const double boundl, const ae
_int_t boundrtype, const double boundr, const real_1d_array &x2, const ae_i
nt_t n2, real_1d_array &y2, real_1d_array &d2, real_1d_array &dd2); | |
| void spline1dconvdiff2cubic(const real_1d_array &x, const real_1d_array &y,
const real_1d_array &x2, real_1d_array &y2, real_1d_array &d2, real_1d_arr
ay &dd2); | | void spline1dconvdiff2cubic(const real_1d_array &x, const real_1d_array &y,
const real_1d_array &x2, real_1d_array &y2, real_1d_array &d2, real_1d_arr
ay &dd2); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
| This subroutine builds Catmull-Rom spline interpolant. | | This subroutine builds Catmull-Rom spline interpolant. | |
| | | | |
|
| INPUT PARAMETERS: | | INPUT PARAMETERS: | |
| X - spline nodes, array[0..N-1]. | | X - spline nodes, array[0..N-1]. | |
| Y - function values, array[0..N-1]. | | Y - function values, array[0..N-1]. | |
| | | | |
| OPTIONAL PARAMETERS: | | | |
| N - points count: | | | |
| * N>=2 | | | |
| * if given, only first N points are used to build splin | | | |
| e | | | |
| * if not given, automatically detected from X/Y sizes | | | |
| (len(X) must be equal to len(Y)) | | | |
| BoundType - boundary condition type: | | | |
| * -1 for periodic boundary condition | | | |
| * 0 for parabolically terminated spline (default) | | | |
| Tension - tension parameter: | | | |
| * tension=0 corresponds to classic Catmull-Rom spline | | | |
| (default) | | | |
| * 0<tension<1 corresponds to more general form - cardin | | | |
| al spline | | | |
| | | | |
| OUTPUT PARAMETERS: | | | |
| C - spline interpolant | | | |
| | | | |
| ORDER OF POINTS | | | |
| | | | |
| Subroutine automatically sorts points, so caller may pass unsorted array. | | | |
| | | | |
| PROBLEMS WITH PERIODIC BOUNDARY CONDITIONS: | | | |
| | | | |
| Problems with periodic boundary conditions have Y[first_point]=Y[last_point | | | |
| ]. | | | |
| However, this subroutine doesn't require you to specify equal values for | | | |
| the first and last points - it automatically forces them to be equal by | | | |
| copying Y[first_point] (corresponds to the leftmost, minimal X[]) to | | | |
| Y[last_point]. However it is recommended to pass consistent values of Y[], | | | |
| i.e. to make Y[first_point]=Y[last_point]. | | | |
| | | | |
| -- ALGLIB PROJECT -- | | | |
| Copyright 23.06.2007 by Bochkanov Sergey | | | |
| *************************************************************************/ | | | |
| void spline1dbuildcatmullrom(const real_1d_array &x, const real_1d_array &y | | | |
| , const ae_int_t n, const ae_int_t boundtype, const double tension, spline1 | | | |
| dinterpolant &c); | | | |
| void spline1dbuildcatmullrom(const real_1d_array &x, const real_1d_array &y | | | |
| , spline1dinterpolant &c); | | | |
| | | | |
| /************************************************************************* | | | |
| This subroutine builds Hermite spline interpolant. | | | |
| | | | |
| INPUT PARAMETERS: | | | |
| X - spline nodes, array[0..N-1] | | | |
| Y - function values, array[0..N-1] | | | |
| D - derivatives, array[0..N-1] | | | |
| N - points count (optional): | | | |
| * N>=2 | | | |
| * if given, only first N points are used to build splin | | | |
| e | | | |
| * if not given, automatically detected from X/Y sizes | | | |
| (len(X) must be equal to len(Y)) | | | |
| | | | |
| OUTPUT PARAMETERS: | | | |
| C - spline interpolant. | | | |
| | | | |
| ORDER OF POINTS | | | |
| | | | |
| Subroutine automatically sorts points, so caller may pass unsorted array. | | | |
| | | | |
| -- ALGLIB PROJECT -- | | | |
| Copyright 23.06.2007 by Bochkanov Sergey | | | |
| *************************************************************************/ | | | |
| void spline1dbuildhermite(const real_1d_array &x, const real_1d_array &y, c | | | |
| onst real_1d_array &d, const ae_int_t n, spline1dinterpolant &c); | | | |
| void spline1dbuildhermite(const real_1d_array &x, const real_1d_array &y, c | | | |
| onst real_1d_array &d, spline1dinterpolant &c); | | | |
| | | | |
| /************************************************************************* | | | |
| This subroutine builds Akima spline interpolant | | | |
| | | | |
| INPUT PARAMETERS: | | | |
| X - spline nodes, array[0..N-1] | | | |
| Y - function values, array[0..N-1] | | | |
| N - points count (optional): | | | |
| * N>=5 | | | |
| * if given, only first N points are used to build splin | | | |
| e | | | |
| * if not given, automatically detected from X/Y sizes | | | |
| (len(X) must be equal to len(Y)) | | | |
| | | | |
| OUTPUT PARAMETERS: | | | |
| C - spline interpolant | | | |
| | | | |
| ORDER OF POINTS | | | |
| | | | |
| Subroutine automatically sorts points, so caller may pass unsorted array. | | | |
| | | | |
| -- ALGLIB PROJECT -- | | | |
| Copyright 24.06.2007 by Bochkanov Sergey | | | |
| *************************************************************************/ | | | |
| void spline1dbuildakima(const real_1d_array &x, const real_1d_array &y, con | | | |
| st ae_int_t n, spline1dinterpolant &c); | | | |
| void spline1dbuildakima(const real_1d_array &x, const real_1d_array &y, spl | | | |
| ine1dinterpolant &c); | | | |
| | | | |
| /************************************************************************* | | | |
| Weighted fitting by cubic spline, with constraints on function values or | | | |
| derivatives. | | | |
| | | | |
| Equidistant grid with M-2 nodes on [min(x,xc),max(x,xc)] is used to build | | | |
| basis functions. Basis functions are cubic splines with continuous second | | | |
| derivatives and non-fixed first derivatives at interval ends. Small | | | |
| regularizing term is used when solving constrained tasks (to improve | | | |
| stability). | | | |
| | | | |
| Task is linear, so linear least squares solver is used. Complexity of this | | | |
| computational scheme is O(N*M^2), mostly dominated by least squares solver | | | |
| | | | |
| SEE ALSO | | | |
| Spline1DFitHermiteWC() - fitting by Hermite splines (more flexible, | | | |
| less smooth) | | | |
| Spline1DFitCubic() - "lightweight" fitting by cubic splines, | | | |
| without invididual weights and constraints | | | |
| | | | |
| INPUT PARAMETERS: | | | |
| X - points, array[0..N-1]. | | | |
| Y - function values, array[0..N-1]. | | | |
| W - weights, array[0..N-1] | | | |
| Each summand in square sum of approximation deviations from | | | |
| given values is multiplied by the square of corresponding | | | |
| weight. Fill it by 1's if you don't want to solve weighted | | | |
| task. | | | |
| N - number of points (optional): | | | |
| * N>0 | | | |
| * if given, only first N elements of X/Y/W are processed | | | |
| * if not given, automatically determined from X/Y/W sizes | | | |
| XC - points where spline values/derivatives are constrained, | | | |
| array[0..K-1]. | | | |
| YC - values of constraints, array[0..K-1] | | | |
| DC - array[0..K-1], types of constraints: | | | |
| * DC[i]=0 means that S(XC[i])=YC[i] | | | |
| * DC[i]=1 means that S'(XC[i])=YC[i] | | | |
| SEE BELOW FOR IMPORTANT INFORMATION ON CONSTRAINTS | | | |
| K - number of constraints (optional): | | | |
| * 0<=K<M. | | | |
| * K=0 means no constraints (XC/YC/DC are not used) | | | |
| * if given, only first K elements of XC/YC/DC are used | | | |
| * if not given, automatically determined from XC/YC/DC | | | |
| M - number of basis functions ( = number_of_nodes+2), M>=4. | | | |
| | | | |
| OUTPUT PARAMETERS: | | | |
| Info- same format as in LSFitLinearWC() subroutine. | | | |
| * Info>0 task is solved | | | |
| * Info<=0 an error occured: | | | |
| -4 means inconvergence of internal SVD | | | |
| -3 means inconsistent constraints | | | |
| S - spline interpolant. | | | |
| Rep - report, same format as in LSFitLinearWC() subroutine. | | | |
| Following fields are set: | | | |
| * RMSError rms error on the (X,Y). | | | |
| * AvgError average error on the (X,Y). | | | |
| * AvgRelError average relative error on the non-zero Y | | | |
| * MaxError maximum error | | | |
| NON-WEIGHTED ERRORS ARE CALCULATED | | | |
| | | | |
| IMPORTANT: | | | |
| this subroitine doesn't calculate task's condition number for K<>0. | | | |
| | | | |
| ORDER OF POINTS | | | |
| | | | |
| Subroutine automatically sorts points, so caller may pass unsorted array. | | | |
| | | | |
| SETTING CONSTRAINTS - DANGERS AND OPPORTUNITIES: | | | |
| | | | |
| Setting constraints can lead to undesired results, like ill-conditioned | | | |
| behavior, or inconsistency being detected. From the other side, it allows | | | |
| us to improve quality of the fit. Here we summarize our experience with | | | |
| constrained regression splines: | | | |
| * excessive constraints can be inconsistent. Splines are piecewise cubic | | | |
| functions, and it is easy to create an example, where large number of | | | |
| constraints concentrated in small area will result in inconsistency. | | | |
| Just because spline is not flexible enough to satisfy all of them. And | | | |
| same constraints spread across the [min(x),max(x)] will be perfectly | | | |
| consistent. | | | |
| * the more evenly constraints are spread across [min(x),max(x)], the more | | | |
| chances that they will be consistent | | | |
| * the greater is M (given fixed constraints), the more chances that | | | |
| constraints will be consistent | | | |
| * in the general case, consistency of constraints IS NOT GUARANTEED. | | | |
| * in the several special cases, however, we CAN guarantee consistency. | | | |
| * one of this cases is constraints on the function values AND/OR its | | | |
| derivatives at the interval boundaries. | | | |
| * another special case is ONE constraint on the function value (OR, but | | | |
| not AND, derivative) anywhere in the interval | | | |
| | | | |
| Our final recommendation is to use constraints WHEN AND ONLY WHEN you | | | |
| can't solve your task without them. Anything beyond special cases given | | | |
| above is not guaranteed and may result in inconsistency. | | | |
| | | | |
| -- ALGLIB PROJECT -- | | | |
| Copyright 18.08.2009 by Bochkanov Sergey | | | |
| *************************************************************************/ | | | |
| void spline1dfitcubicwc(const real_1d_array &x, const real_1d_array &y, con | | | |
| st real_1d_array &w, const ae_int_t n, const real_1d_array &xc, const real_ | | | |
| 1d_array &yc, const integer_1d_array &dc, const ae_int_t k, const ae_int_t | | | |
| m, ae_int_t &info, spline1dinterpolant &s, spline1dfitreport &rep); | | | |
| void spline1dfitcubicwc(const real_1d_array &x, const real_1d_array &y, con | | | |
| st real_1d_array &w, const real_1d_array &xc, const real_1d_array &yc, cons | | | |
| t integer_1d_array &dc, const ae_int_t m, ae_int_t &info, spline1dinterpola | | | |
| nt &s, spline1dfitreport &rep); | | | |
| | | | |
| /************************************************************************* | | | |
| Weighted fitting by Hermite spline, with constraints on function values | | | |
| or first derivatives. | | | |
| | | | |
| Equidistant grid with M nodes on [min(x,xc),max(x,xc)] is used to build | | | |
| basis functions. Basis functions are Hermite splines. Small regularizing | | | |
| term is used when solving constrained tasks (to improve stability). | | | |
| | | | |
| Task is linear, so linear least squares solver is used. Complexity of this | | | |
| computational scheme is O(N*M^2), mostly dominated by least squares solver | | | |
| | | | |
| SEE ALSO | | | |
| Spline1DFitCubicWC() - fitting by Cubic splines (less flexible, | | | |
| more smooth) | | | |
| Spline1DFitHermite() - "lightweight" Hermite fitting, without | | | |
| invididual weights and constraints | | | |
| | | | |
| INPUT PARAMETERS: | | | |
| X - points, array[0..N-1]. | | | |
| Y - function values, array[0..N-1]. | | | |
| W - weights, array[0..N-1] | | | |
| Each summand in square sum of approximation deviations from | | | |
| given values is multiplied by the square of corresponding | | | |
| weight. Fill it by 1's if you don't want to solve weighted | | | |
| task. | | | |
| N - number of points (optional): | | | |
| * N>0 | | | |
| * if given, only first N elements of X/Y/W are processed | | | |
| * if not given, automatically determined from X/Y/W sizes | | | |
| XC - points where spline values/derivatives are constrained, | | | |
| array[0..K-1]. | | | |
| YC - values of constraints, array[0..K-1] | | | |
| DC - array[0..K-1], types of constraints: | | | |
| * DC[i]=0 means that S(XC[i])=YC[i] | | | |
| * DC[i]=1 means that S'(XC[i])=YC[i] | | | |
| SEE BELOW FOR IMPORTANT INFORMATION ON CONSTRAINTS | | | |
| K - number of constraints (optional): | | | |
| * 0<=K<M. | | | |
| * K=0 means no constraints (XC/YC/DC are not used) | | | |
| * if given, only first K elements of XC/YC/DC are used | | | |
| * if not given, automatically determined from XC/YC/DC | | | |
| M - number of basis functions (= 2 * number of nodes), | | | |
| M>=4, | | | |
| M IS EVEN! | | | |
| | | | |
| OUTPUT PARAMETERS: | | | |
| Info- same format as in LSFitLinearW() subroutine: | | | |
| * Info>0 task is solved | | | |
| * Info<=0 an error occured: | | | |
| -4 means inconvergence of internal SVD | | | |
| -3 means inconsistent constraints | | | |
| -2 means odd M was passed (which is not supported) | | | |
| -1 means another errors in parameters passed | | | |
| (N<=0, for example) | | | |
| S - spline interpolant. | | | |
| Rep - report, same format as in LSFitLinearW() subroutine. | | | |
| Following fields are set: | | | |
| * RMSError rms error on the (X,Y). | | | |
| * AvgError average error on the (X,Y). | | | |
| * AvgRelError average relative error on the non-zero Y | | | |
| * MaxError maximum error | | | |
| NON-WEIGHTED ERRORS ARE CALCULATED | | | |
| | | | |
| IMPORTANT: | | | |
| this subroitine doesn't calculate task's condition number for K<>0. | | | |
| | | | |
|
| IMPORTANT: | | OPTIONAL PARAMETERS: | |
| this subroitine supports only even M's | | N - points count: | |
| | | * N>=2 | |
| | | * if given, only first N points are used to build splin | |
| | | e | |
| | | * if not given, automatically detected from X/Y sizes | |
| | | (len(X) must be equal to len(Y)) | |
| | | BoundType - boundary condition type: | |
| | | * -1 for periodic boundary condition | |
| | | * 0 for parabolically terminated spline (default) | |
| | | Tension - tension parameter: | |
| | | * tension=0 corresponds to classic Catmull-Rom spline | |
| | | (default) | |
| | | * 0<tension<1 corresponds to more general form - cardin | |
| | | al spline | |
| | | | |
| | | OUTPUT PARAMETERS: | |
| | | C - spline interpolant | |
| | | | |
| ORDER OF POINTS | | ORDER OF POINTS | |
| | | | |
| Subroutine automatically sorts points, so caller may pass unsorted array. | | Subroutine automatically sorts points, so caller may pass unsorted array. | |
| | | | |
|
| SETTING CONSTRAINTS - DANGERS AND OPPORTUNITIES: | | PROBLEMS WITH PERIODIC BOUNDARY CONDITIONS: | |
| | | | |
| Setting constraints can lead to undesired results, like ill-conditioned | | | |
| behavior, or inconsistency being detected. From the other side, it allows | | | |
| us to improve quality of the fit. Here we summarize our experience with | | | |
| constrained regression splines: | | | |
| * excessive constraints can be inconsistent. Splines are piecewise cubic | | | |
| functions, and it is easy to create an example, where large number of | | | |
| constraints concentrated in small area will result in inconsistency. | | | |
| Just because spline is not flexible enough to satisfy all of them. And | | | |
| same constraints spread across the [min(x),max(x)] will be perfectly | | | |
| consistent. | | | |
| * the more evenly constraints are spread across [min(x),max(x)], the more | | | |
| chances that they will be consistent | | | |
| * the greater is M (given fixed constraints), the more chances that | | | |
| constraints will be consistent | | | |
| * in the general case, consistency of constraints is NOT GUARANTEED. | | | |
| * in the several special cases, however, we can guarantee consistency. | | | |
| * one of this cases is M>=4 and constraints on the function value | | | |
| (AND/OR its derivative) at the interval boundaries. | | | |
| * another special case is M>=4 and ONE constraint on the function value | | | |
| (OR, BUT NOT AND, derivative) anywhere in [min(x),max(x)] | | | |
| | | | |
|
| Our final recommendation is to use constraints WHEN AND ONLY when you | | Problems with periodic boundary conditions have Y[first_point]=Y[last_point | |
| can't solve your task without them. Anything beyond special cases given | | ]. | |
| above is not guaranteed and may result in inconsistency. | | However, this subroutine doesn't require you to specify equal values for | |
| | | the first and last points - it automatically forces them to be equal by | |
| | | copying Y[first_point] (corresponds to the leftmost, minimal X[]) to | |
| | | Y[last_point]. However it is recommended to pass consistent values of Y[], | |
| | | i.e. to make Y[first_point]=Y[last_point]. | |
| | | | |
| -- ALGLIB PROJECT -- | | -- ALGLIB PROJECT -- | |
|
| Copyright 18.08.2009 by Bochkanov Sergey | | Copyright 23.06.2007 by Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
|
| void spline1dfithermitewc(const real_1d_array &x, const real_1d_array &y, c | | void spline1dbuildcatmullrom(const real_1d_array &x, const real_1d_array &y | |
| onst real_1d_array &w, const ae_int_t n, const real_1d_array &xc, const rea | | , const ae_int_t n, const ae_int_t boundtype, const double tension, spline1 | |
| l_1d_array &yc, const integer_1d_array &dc, const ae_int_t k, const ae_int_ | | dinterpolant &c); | |
| t m, ae_int_t &info, spline1dinterpolant &s, spline1dfitreport &rep); | | void spline1dbuildcatmullrom(const real_1d_array &x, const real_1d_array &y | |
| void spline1dfithermitewc(const real_1d_array &x, const real_1d_array &y, c | | , spline1dinterpolant &c); | |
| onst real_1d_array &w, const real_1d_array &xc, const real_1d_array &yc, co | | | |
| nst integer_1d_array &dc, const ae_int_t m, ae_int_t &info, spline1dinterpo | | | |
| lant &s, spline1dfitreport &rep); | | | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
|
| Least squares fitting by cubic spline. | | This subroutine builds Hermite spline interpolant. | |
| | | | |
|
| This subroutine is "lightweight" alternative for more complex and feature- | | INPUT PARAMETERS: | |
| rich Spline1DFitCubicWC(). See Spline1DFitCubicWC() for more information | | X - spline nodes, array[0..N-1] | |
| about subroutine parameters (we don't duplicate it here because of length) | | Y - function values, array[0..N-1] | |
| | | D - derivatives, array[0..N-1] | |
| | | N - points count (optional): | |
| | | * N>=2 | |
| | | * if given, only first N points are used to build splin | |
| | | e | |
| | | * if not given, automatically detected from X/Y sizes | |
| | | (len(X) must be equal to len(Y)) | |
| | | | |
| | | OUTPUT PARAMETERS: | |
| | | C - spline interpolant. | |
| | | | |
| | | ORDER OF POINTS | |
| | | | |
| | | Subroutine automatically sorts points, so caller may pass unsorted array. | |
| | | | |
| -- ALGLIB PROJECT -- | | -- ALGLIB PROJECT -- | |
|
| Copyright 18.08.2009 by Bochkanov Sergey | | Copyright 23.06.2007 by Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
|
| void spline1dfitcubic(const real_1d_array &x, const real_1d_array &y, const | | void spline1dbuildhermite(const real_1d_array &x, const real_1d_array &y, c | |
| ae_int_t n, const ae_int_t m, ae_int_t &info, spline1dinterpolant &s, spli | | onst real_1d_array &d, const ae_int_t n, spline1dinterpolant &c); | |
| ne1dfitreport &rep); | | void spline1dbuildhermite(const real_1d_array &x, const real_1d_array &y, c | |
| void spline1dfitcubic(const real_1d_array &x, const real_1d_array &y, const | | onst real_1d_array &d, spline1dinterpolant &c); | |
| ae_int_t m, ae_int_t &info, spline1dinterpolant &s, spline1dfitreport &rep | | | |
| ); | | | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
|
| Least squares fitting by Hermite spline. | | This subroutine builds Akima spline interpolant | |
| | | | |
|
| This subroutine is "lightweight" alternative for more complex and feature- | | INPUT PARAMETERS: | |
| rich Spline1DFitHermiteWC(). See Spline1DFitHermiteWC() description for | | X - spline nodes, array[0..N-1] | |
| more information about subroutine parameters (we don't duplicate it here | | Y - function values, array[0..N-1] | |
| because of length). | | N - points count (optional): | |
| | | * N>=5 | |
| | | * if given, only first N points are used to build splin | |
| | | e | |
| | | * if not given, automatically detected from X/Y sizes | |
| | | (len(X) must be equal to len(Y)) | |
| | | | |
| | | OUTPUT PARAMETERS: | |
| | | C - spline interpolant | |
| | | | |
| | | ORDER OF POINTS | |
| | | | |
| | | Subroutine automatically sorts points, so caller may pass unsorted array. | |
| | | | |
| -- ALGLIB PROJECT -- | | -- ALGLIB PROJECT -- | |
|
| Copyright 18.08.2009 by Bochkanov Sergey | | Copyright 24.06.2007 by Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
|
| void spline1dfithermite(const real_1d_array &x, const real_1d_array &y, con | | void spline1dbuildakima(const real_1d_array &x, const real_1d_array &y, con | |
| st ae_int_t n, const ae_int_t m, ae_int_t &info, spline1dinterpolant &s, sp | | st ae_int_t n, spline1dinterpolant &c); | |
| line1dfitreport &rep); | | void spline1dbuildakima(const real_1d_array &x, const real_1d_array &y, spl | |
| void spline1dfithermite(const real_1d_array &x, const real_1d_array &y, con | | ine1dinterpolant &c); | |
| st ae_int_t m, ae_int_t &info, spline1dinterpolant &s, spline1dfitreport &r | | | |
| ep); | | | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
| This subroutine calculates the value of the spline at the given point X. | | This subroutine calculates the value of the spline at the given point X. | |
| | | | |
| INPUT PARAMETERS: | | INPUT PARAMETERS: | |
| C - spline interpolant | | C - spline interpolant | |
| X - point | | X - point | |
| | | | |
| Result: | | Result: | |
| S(x) | | S(x) | |
| | | | |
| skipping to change at line 2437 | | skipping to change at line 2700 | |
|---|
| Tbl[I,3] = C1 | | Tbl[I,3] = C1 | |
| Tbl[I,4] = C2 | | Tbl[I,4] = C2 | |
| Tbl[I,5] = C3 | | Tbl[I,5] = C3 | |
| On [x[i], x[i+1]] spline is equals to: | | On [x[i], x[i+1]] spline is equals to: | |
| S(x) = C0 + C1*t + C2*t^2 + C3*t^3 | | S(x) = C0 + C1*t + C2*t^2 + C3*t^3 | |
| t = x-x[i] | | t = x-x[i] | |
| | | | |
| -- ALGLIB PROJECT -- | | -- ALGLIB PROJECT -- | |
| Copyright 29.06.2007 by Bochkanov Sergey | | Copyright 29.06.2007 by Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
|
| void spline1dunpack(const spline1dinterpolant &c, ae_int_t &n, real_2d_arra | | void spline1dunpack(const spline1dinterpolant &c, ae_int_t &n, real_2d_arra | |
| y &tbl); | | y &tbl); | |
| | | | |
| | | /************************************************************************* | |
| | | This subroutine performs linear transformation of the spline argument. | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | C - spline interpolant. | |
| | | A, B- transformation coefficients: x = A*t + B | |
| | | Result: | |
| | | C - transformed spline | |
| | | | |
| | | -- ALGLIB PROJECT -- | |
| | | Copyright 30.06.2007 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void spline1dlintransx(const spline1dinterpolant &c, const double a, const | |
| | | double b); | |
| | | | |
| | | /************************************************************************* | |
| | | This subroutine performs linear transformation of the spline. | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | C - spline interpolant. | |
| | | A, B- transformation coefficients: S2(x) = A*S(x) + B | |
| | | Result: | |
| | | C - transformed spline | |
| | | | |
| | | -- ALGLIB PROJECT -- | |
| | | Copyright 30.06.2007 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void spline1dlintransy(const spline1dinterpolant &c, const double a, const | |
| | | double b); | |
| | | | |
| | | /************************************************************************* | |
| | | This subroutine integrates the spline. | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | C - spline interpolant. | |
| | | X - right bound of the integration interval [a, x], | |
| | | here 'a' denotes min(x[]) | |
| | | Result: | |
| | | integral(S(t)dt,a,x) | |
| | | | |
| | | -- ALGLIB PROJECT -- | |
| | | Copyright 23.06.2007 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | double spline1dintegrate(const spline1dinterpolant &c, const double x); | |
| | | | |
| | | /************************************************************************* | |
| | | This subroutine builds bilinear spline coefficients table. | |
| | | | |
| | | Input parameters: | |
| | | X - spline abscissas, array[0..N-1] | |
| | | Y - spline ordinates, array[0..M-1] | |
| | | F - function values, array[0..M-1,0..N-1] | |
| | | M,N - grid size, M>=2, N>=2 | |
| | | | |
| | | Output parameters: | |
| | | C - spline interpolant | |
| | | | |
| | | -- ALGLIB PROJECT -- | |
| | | Copyright 05.07.2007 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void spline2dbuildbilinear(const real_1d_array &x, const real_1d_array &y, | |
| | | const real_2d_array &f, const ae_int_t m, const ae_int_t n, spline2dinterpo | |
| | | lant &c); | |
| | | | |
| | | /************************************************************************* | |
| | | This subroutine builds bicubic spline coefficients table. | |
| | | | |
| | | Input parameters: | |
| | | X - spline abscissas, array[0..N-1] | |
| | | Y - spline ordinates, array[0..M-1] | |
| | | F - function values, array[0..M-1,0..N-1] | |
| | | M,N - grid size, M>=2, N>=2 | |
| | | | |
| | | Output parameters: | |
| | | C - spline interpolant | |
| | | | |
| | | -- ALGLIB PROJECT -- | |
| | | Copyright 05.07.2007 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void spline2dbuildbicubic(const real_1d_array &x, const real_1d_array &y, c | |
| | | onst real_2d_array &f, const ae_int_t m, const ae_int_t n, spline2dinterpol | |
| | | ant &c); | |
| | | | |
| | | /************************************************************************* | |
| | | This subroutine calculates the value of the bilinear or bicubic spline at | |
| | | the given point X. | |
| | | | |
| | | Input parameters: | |
| | | C - coefficients table. | |
| | | Built by BuildBilinearSpline or BuildBicubicSpline. | |
| | | X, Y- point | |
| | | | |
| | | Result: | |
| | | S(x,y) | |
| | | | |
| | | -- ALGLIB PROJECT -- | |
| | | Copyright 05.07.2007 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | double spline2dcalc(const spline2dinterpolant &c, const double x, const dou | |
| | | ble y); | |
| | | | |
| | | /************************************************************************* | |
| | | This subroutine calculates the value of the bilinear or bicubic spline at | |
| | | the given point X and its derivatives. | |
| | | | |
| | | Input parameters: | |
| | | C - spline interpolant. | |
| | | X, Y- point | |
| | | | |
| | | Output parameters: | |
| | | F - S(x,y) | |
| | | FX - dS(x,y)/dX | |
| | | FY - dS(x,y)/dY | |
| | | FXY - d2S(x,y)/dXdY | |
| | | | |
| | | -- ALGLIB PROJECT -- | |
| | | Copyright 05.07.2007 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void spline2ddiff(const spline2dinterpolant &c, const double x, const doubl | |
| | | e y, double &f, double &fx, double &fy, double &fxy); | |
| | | | |
| | | /************************************************************************* | |
| | | This subroutine unpacks two-dimensional spline into the coefficients table | |
| | | | |
| | | Input parameters: | |
| | | C - spline interpolant. | |
| | | | |
| | | Result: | |
| | | M, N- grid size (x-axis and y-axis) | |
| | | Tbl - coefficients table, unpacked format, | |
| | | [0..(N-1)*(M-1)-1, 0..19]. | |
| | | For I = 0...M-2, J=0..N-2: | |
| | | K = I*(N-1)+J | |
| | | Tbl[K,0] = X[j] | |
| | | Tbl[K,1] = X[j+1] | |
| | | Tbl[K,2] = Y[i] | |
| | | Tbl[K,3] = Y[i+1] | |
| | | Tbl[K,4] = C00 | |
| | | Tbl[K,5] = C01 | |
| | | Tbl[K,6] = C02 | |
| | | Tbl[K,7] = C03 | |
| | | Tbl[K,8] = C10 | |
| | | Tbl[K,9] = C11 | |
| | | ... | |
| | | Tbl[K,19] = C33 | |
| | | On each grid square spline is equals to: | |
| | | S(x) = SUM(c[i,j]*(x^i)*(y^j), i=0..3, j=0..3) | |
| | | t = x-x[j] | |
| | | u = y-y[i] | |
| | | | |
| | | -- ALGLIB PROJECT -- | |
| | | Copyright 29.06.2007 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void spline2dunpack(const spline2dinterpolant &c, ae_int_t &m, ae_int_t &n, | |
| | | real_2d_array &tbl); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
| This subroutine performs linear transformation of the spline argument. | | This subroutine performs linear transformation of the spline argument. | |
| | | | |
|
| INPUT PARAMETERS: | | Input parameters: | |
| C - spline interpolant. | | C - spline interpolant | |
| A, B- transformation coefficients: x = A*t + B | | AX, BX - transformation coefficients: x = A*t + B | |
| | | AY, BY - transformation coefficients: y = A*u + B | |
| Result: | | Result: | |
| C - transformed spline | | C - transformed spline | |
| | | | |
| -- ALGLIB PROJECT -- | | -- ALGLIB PROJECT -- | |
| Copyright 30.06.2007 by Bochkanov Sergey | | Copyright 30.06.2007 by Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
|
| void spline1dlintransx(const spline1dinterpolant &c, const double a, const
double b); | | void spline2dlintransxy(const spline2dinterpolant &c, const double ax, cons
t double bx, const double ay, const double by); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
| This subroutine performs linear transformation of the spline. | | This subroutine performs linear transformation of the spline. | |
| | | | |
|
| INPUT PARAMETERS: | | Input parameters: | |
| C - spline interpolant. | | C - spline interpolant. | |
|
| A, B- transformation coefficients: S2(x) = A*S(x) + B | | A, B- transformation coefficients: S2(x,y) = A*S(x,y) + B | |
| Result: | | | |
| | | Output parameters: | |
| C - transformed spline | | C - transformed spline | |
| | | | |
| -- ALGLIB PROJECT -- | | -- ALGLIB PROJECT -- | |
| Copyright 30.06.2007 by Bochkanov Sergey | | Copyright 30.06.2007 by Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
|
| void spline1dlintransy(const spline1dinterpolant &c, const double a, const
double b); | | void spline2dlintransf(const spline2dinterpolant &c, const double a, const
double b); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
|
| This subroutine integrates the spline. | | Bicubic spline resampling | |
| | | | |
|
| INPUT PARAMETERS: | | Input parameters: | |
| C - spline interpolant. | | A - function values at the old grid, | |
| X - right bound of the integration interval [a, x], | | array[0..OldHeight-1, 0..OldWidth-1] | |
| here 'a' denotes min(x[]) | | OldHeight - old grid height, OldHeight>1 | |
| Result: | | OldWidth - old grid width, OldWidth>1 | |
| integral(S(t)dt,a,x) | | NewHeight - new grid height, NewHeight>1 | |
| | | NewWidth - new grid width, NewWidth>1 | |
| | | | |
|
| -- ALGLIB PROJECT -- | | Output parameters: | |
| Copyright 23.06.2007 by Bochkanov Sergey | | B - function values at the new grid, | |
| | | array[0..NewHeight-1, 0..NewWidth-1] | |
| | | | |
| | | -- ALGLIB routine -- | |
| | | 15 May, 2007 | |
| | | Copyright by Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
|
| double spline1dintegrate(const spline1dinterpolant &c, const double x); | | void spline2dresamplebicubic(const real_2d_array &a, const ae_int_t oldheig | |
| | | ht, const ae_int_t oldwidth, real_2d_array &b, const ae_int_t newheight, co | |
| | | nst ae_int_t newwidth); | |
| | | | |
| | | /************************************************************************* | |
| | | Bilinear spline resampling | |
| | | | |
| | | Input parameters: | |
| | | A - function values at the old grid, | |
| | | array[0..OldHeight-1, 0..OldWidth-1] | |
| | | OldHeight - old grid height, OldHeight>1 | |
| | | OldWidth - old grid width, OldWidth>1 | |
| | | NewHeight - new grid height, NewHeight>1 | |
| | | NewWidth - new grid width, NewWidth>1 | |
| | | | |
| | | Output parameters: | |
| | | B - function values at the new grid, | |
| | | array[0..NewHeight-1, 0..NewWidth-1] | |
| | | | |
| | | -- ALGLIB routine -- | |
| | | 09.07.2007 | |
| | | Copyright by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void spline2dresamplebilinear(const real_2d_array &a, const ae_int_t oldhei | |
| | | ght, const ae_int_t oldwidth, real_2d_array &b, const ae_int_t newheight, c | |
| | | onst ae_int_t newwidth); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
| IDW interpolation | | IDW interpolation | |
| | | | |
| INPUT PARAMETERS: | | INPUT PARAMETERS: | |
| Z - IDW interpolant built with one of model building | | Z - IDW interpolant built with one of model building | |
| subroutines. | | subroutines. | |
| X - array[0..NX-1], interpolation point | | X - array[0..NX-1], interpolation point | |
| | | | |
| Result: | | Result: | |
| | | | |
| skipping to change at line 3015 | | skipping to change at line 3454 | |
|---|
| 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 SECTION CONTAINS COMPUTATIONAL CORE DECLARATIONS (FUNCTIONS) | | // THIS SECTION CONTAINS COMPUTATIONAL CORE DECLARATIONS (FUNCTIONS) | |
| // | | // | |
| ///////////////////////////////////////////////////////////////////////// | | ///////////////////////////////////////////////////////////////////////// | |
| namespace alglib_impl | | namespace alglib_impl | |
| { | | { | |
|
| | | void polynomialfit(/* Real */ ae_vector* x, | |
| | | /* Real */ ae_vector* y, | |
| | | ae_int_t n, | |
| | | ae_int_t m, | |
| | | ae_int_t* info, | |
| | | barycentricinterpolant* p, | |
| | | polynomialfitreport* rep, | |
| | | ae_state *_state); | |
| | | void polynomialfitwc(/* Real */ ae_vector* x, | |
| | | /* Real */ ae_vector* y, | |
| | | /* Real */ ae_vector* w, | |
| | | ae_int_t n, | |
| | | /* Real */ ae_vector* xc, | |
| | | /* Real */ ae_vector* yc, | |
| | | /* Integer */ ae_vector* dc, | |
| | | ae_int_t k, | |
| | | ae_int_t m, | |
| | | ae_int_t* info, | |
| | | barycentricinterpolant* p, | |
| | | polynomialfitreport* rep, | |
| | | ae_state *_state); | |
| | | void barycentricfitfloaterhormannwc(/* Real */ ae_vector* x, | |
| | | /* Real */ ae_vector* y, | |
| | | /* Real */ ae_vector* w, | |
| | | ae_int_t n, | |
| | | /* Real */ ae_vector* xc, | |
| | | /* Real */ ae_vector* yc, | |
| | | /* Integer */ ae_vector* dc, | |
| | | ae_int_t k, | |
| | | ae_int_t m, | |
| | | ae_int_t* info, | |
| | | barycentricinterpolant* b, | |
| | | barycentricfitreport* rep, | |
| | | ae_state *_state); | |
| | | void barycentricfitfloaterhormann(/* Real */ ae_vector* x, | |
| | | /* Real */ ae_vector* y, | |
| | | ae_int_t n, | |
| | | ae_int_t m, | |
| | | ae_int_t* info, | |
| | | barycentricinterpolant* b, | |
| | | barycentricfitreport* rep, | |
| | | ae_state *_state); | |
| | | void spline1dfitpenalized(/* Real */ ae_vector* x, | |
| | | /* Real */ ae_vector* y, | |
| | | ae_int_t n, | |
| | | ae_int_t m, | |
| | | double rho, | |
| | | ae_int_t* info, | |
| | | spline1dinterpolant* s, | |
| | | spline1dfitreport* rep, | |
| | | ae_state *_state); | |
| | | void spline1dfitpenalizedw(/* Real */ ae_vector* x, | |
| | | /* Real */ ae_vector* y, | |
| | | /* Real */ ae_vector* w, | |
| | | ae_int_t n, | |
| | | ae_int_t m, | |
| | | double rho, | |
| | | ae_int_t* info, | |
| | | spline1dinterpolant* s, | |
| | | spline1dfitreport* rep, | |
| | | ae_state *_state); | |
| | | void spline1dfitcubicwc(/* Real */ ae_vector* x, | |
| | | /* Real */ ae_vector* y, | |
| | | /* Real */ ae_vector* w, | |
| | | ae_int_t n, | |
| | | /* Real */ ae_vector* xc, | |
| | | /* Real */ ae_vector* yc, | |
| | | /* Integer */ ae_vector* dc, | |
| | | ae_int_t k, | |
| | | ae_int_t m, | |
| | | ae_int_t* info, | |
| | | spline1dinterpolant* s, | |
| | | spline1dfitreport* rep, | |
| | | ae_state *_state); | |
| | | void spline1dfithermitewc(/* Real */ ae_vector* x, | |
| | | /* Real */ ae_vector* y, | |
| | | /* Real */ ae_vector* w, | |
| | | ae_int_t n, | |
| | | /* Real */ ae_vector* xc, | |
| | | /* Real */ ae_vector* yc, | |
| | | /* Integer */ ae_vector* dc, | |
| | | ae_int_t k, | |
| | | ae_int_t m, | |
| | | ae_int_t* info, | |
| | | spline1dinterpolant* s, | |
| | | spline1dfitreport* rep, | |
| | | ae_state *_state); | |
| | | void spline1dfitcubic(/* Real */ ae_vector* x, | |
| | | /* Real */ ae_vector* y, | |
| | | ae_int_t n, | |
| | | ae_int_t m, | |
| | | ae_int_t* info, | |
| | | spline1dinterpolant* s, | |
| | | spline1dfitreport* rep, | |
| | | ae_state *_state); | |
| | | void spline1dfithermite(/* Real */ ae_vector* x, | |
| | | /* Real */ ae_vector* y, | |
| | | ae_int_t n, | |
| | | ae_int_t m, | |
| | | ae_int_t* info, | |
| | | spline1dinterpolant* s, | |
| | | spline1dfitreport* rep, | |
| | | ae_state *_state); | |
| void lsfitlinearw(/* Real */ ae_vector* y, | | void lsfitlinearw(/* Real */ ae_vector* y, | |
| /* Real */ ae_vector* w, | | /* Real */ ae_vector* w, | |
| /* Real */ ae_matrix* fmatrix, | | /* Real */ ae_matrix* fmatrix, | |
| ae_int_t n, | | ae_int_t n, | |
| ae_int_t m, | | ae_int_t m, | |
| ae_int_t* info, | | ae_int_t* info, | |
| /* Real */ ae_vector* c, | | /* Real */ ae_vector* c, | |
| lsfitreport* rep, | | lsfitreport* rep, | |
| ae_state *_state); | | ae_state *_state); | |
| void lsfitlinearwc(/* Real */ ae_vector* y, | | void lsfitlinearwc(/* Real */ ae_vector* y, | |
| | | | |
| skipping to change at line 3053 | | skipping to change at line 3595 | |
|---|
| void lsfitlinearc(/* Real */ ae_vector* y, | | void lsfitlinearc(/* Real */ ae_vector* y, | |
| /* Real */ ae_matrix* fmatrix, | | /* Real */ ae_matrix* fmatrix, | |
| /* Real */ ae_matrix* cmatrix, | | /* Real */ ae_matrix* cmatrix, | |
| ae_int_t n, | | ae_int_t n, | |
| ae_int_t m, | | ae_int_t m, | |
| ae_int_t k, | | ae_int_t k, | |
| ae_int_t* info, | | ae_int_t* info, | |
| /* Real */ ae_vector* c, | | /* Real */ ae_vector* c, | |
| lsfitreport* rep, | | lsfitreport* rep, | |
| ae_state *_state); | | ae_state *_state); | |
|
| | | void lsfitcreatewf(/* Real */ ae_matrix* x, | |
| | | /* Real */ ae_vector* y, | |
| | | /* Real */ ae_vector* w, | |
| | | /* Real */ ae_vector* c, | |
| | | ae_int_t n, | |
| | | ae_int_t m, | |
| | | ae_int_t k, | |
| | | double diffstep, | |
| | | lsfitstate* state, | |
| | | ae_state *_state); | |
| | | void lsfitcreatef(/* Real */ ae_matrix* x, | |
| | | /* Real */ ae_vector* y, | |
| | | /* Real */ ae_vector* c, | |
| | | ae_int_t n, | |
| | | ae_int_t m, | |
| | | ae_int_t k, | |
| | | double diffstep, | |
| | | lsfitstate* state, | |
| | | ae_state *_state); | |
| void lsfitcreatewfg(/* Real */ ae_matrix* x, | | void lsfitcreatewfg(/* Real */ ae_matrix* x, | |
| /* Real */ ae_vector* y, | | /* Real */ ae_vector* y, | |
| /* Real */ ae_vector* w, | | /* Real */ ae_vector* w, | |
| /* Real */ ae_vector* c, | | /* Real */ ae_vector* c, | |
| ae_int_t n, | | ae_int_t n, | |
| ae_int_t m, | | ae_int_t m, | |
| ae_int_t k, | | ae_int_t k, | |
| ae_bool cheapfg, | | ae_bool cheapfg, | |
| lsfitstate* state, | | lsfitstate* state, | |
| ae_state *_state); | | ae_state *_state); | |
| | | | |
| skipping to change at line 3100 | | skipping to change at line 3661 | |
|---|
| double epsf, | | double epsf, | |
| double epsx, | | double epsx, | |
| ae_int_t maxits, | | ae_int_t maxits, | |
| ae_state *_state); | | ae_state *_state); | |
| void lsfitsetstpmax(lsfitstate* state, double stpmax, ae_state *_state); | | void lsfitsetstpmax(lsfitstate* state, double stpmax, ae_state *_state); | |
| void lsfitsetxrep(lsfitstate* state, ae_bool needxrep, ae_state *_state); | | void lsfitsetxrep(lsfitstate* state, ae_bool needxrep, ae_state *_state); | |
| ae_bool lsfititeration(lsfitstate* state, ae_state *_state); | | ae_bool lsfititeration(lsfitstate* state, ae_state *_state); | |
| void lsfitresults(lsfitstate* state, | | void lsfitresults(lsfitstate* state, | |
| ae_int_t* info, | | ae_int_t* info, | |
| /* Real */ ae_vector* c, | | /* Real */ ae_vector* c, | |
|
| lsfitreport* rep, | | lsfitreport* rep, | |
| ae_state *_state); | | | |
| void lsfitscalexy(/* Real */ ae_vector* x, | | | |
| /* Real */ ae_vector* y, | | | |
| ae_int_t n, | | | |
| /* Real */ ae_vector* xc, | | | |
| /* Real */ ae_vector* yc, | | | |
| /* Integer */ ae_vector* dc, | | | |
| ae_int_t k, | | | |
| double* xa, | | | |
| double* xb, | | | |
| double* sa, | | | |
| double* sb, | | | |
| /* Real */ ae_vector* xoriginal, | | | |
| /* Real */ ae_vector* yoriginal, | | | |
| ae_state *_state); | | | |
| ae_bool _lsfitreport_init(lsfitreport* p, ae_state *_state, ae_bool make_au | | | |
| tomatic); | | | |
| ae_bool _lsfitreport_init_copy(lsfitreport* dst, lsfitreport* src, ae_state | | | |
| *_state, ae_bool make_automatic); | | | |
| void _lsfitreport_clear(lsfitreport* p); | | | |
| ae_bool _lsfitstate_init(lsfitstate* p, ae_state *_state, ae_bool make_auto | | | |
| matic); | | | |
| ae_bool _lsfitstate_init_copy(lsfitstate* dst, lsfitstate* src, ae_state *_ | | | |
| state, ae_bool make_automatic); | | | |
| void _lsfitstate_clear(lsfitstate* p); | | | |
| void polynomialbuild(/* Real */ ae_vector* x, | | | |
| /* Real */ ae_vector* y, | | | |
| ae_int_t n, | | | |
| barycentricinterpolant* p, | | | |
| ae_state *_state); | | | |
| void polynomialbuildeqdist(double a, | | | |
| double b, | | | |
| /* Real */ ae_vector* y, | | | |
| ae_int_t n, | | | |
| barycentricinterpolant* p, | | | |
| ae_state *_state); | | | |
| void polynomialbuildcheb1(double a, | | | |
| double b, | | | |
| /* Real */ ae_vector* y, | | | |
| ae_int_t n, | | | |
| barycentricinterpolant* p, | | | |
| ae_state *_state); | | | |
| void polynomialbuildcheb2(double a, | | | |
| double b, | | | |
| /* Real */ ae_vector* y, | | | |
| ae_int_t n, | | | |
| barycentricinterpolant* p, | | | |
| ae_state *_state); | | | |
| double polynomialcalceqdist(double a, | | | |
| double b, | | | |
| /* Real */ ae_vector* f, | | | |
| ae_int_t n, | | | |
| double t, | | | |
| ae_state *_state); | | | |
| double polynomialcalccheb1(double a, | | | |
| double b, | | | |
| /* Real */ ae_vector* f, | | | |
| ae_int_t n, | | | |
| double t, | | | |
| ae_state *_state); | | | |
| double polynomialcalccheb2(double a, | | | |
| double b, | | | |
| /* Real */ ae_vector* f, | | | |
| ae_int_t n, | | | |
| double t, | | | |
| ae_state *_state); | | | |
| void polynomialfit(/* Real */ ae_vector* x, | | | |
| /* Real */ ae_vector* y, | | | |
| ae_int_t n, | | | |
| ae_int_t m, | | | |
| ae_int_t* info, | | | |
| barycentricinterpolant* p, | | | |
| polynomialfitreport* rep, | | | |
| ae_state *_state); | | ae_state *_state); | |
|
| void polynomialfitwc(/* Real */ ae_vector* x, | | void lsfitscalexy(/* Real */ ae_vector* x, | |
| /* Real */ ae_vector* y, | | /* Real */ ae_vector* y, | |
| /* Real */ ae_vector* w, | | /* Real */ ae_vector* w, | |
| ae_int_t n, | | ae_int_t n, | |
| /* Real */ ae_vector* xc, | | /* Real */ ae_vector* xc, | |
| /* Real */ ae_vector* yc, | | /* Real */ ae_vector* yc, | |
| /* Integer */ ae_vector* dc, | | /* Integer */ ae_vector* dc, | |
| ae_int_t k, | | ae_int_t k, | |
|
| ae_int_t m, | | double* xa, | |
| ae_int_t* info, | | double* xb, | |
| barycentricinterpolant* p, | | double* sa, | |
| polynomialfitreport* rep, | | double* sb, | |
| | | /* Real */ ae_vector* xoriginal, | |
| | | /* Real */ ae_vector* yoriginal, | |
| ae_state *_state); | | ae_state *_state); | |
| ae_bool _polynomialfitreport_init(polynomialfitreport* p, ae_state *_state,
ae_bool make_automatic); | | ae_bool _polynomialfitreport_init(polynomialfitreport* p, ae_state *_state,
ae_bool make_automatic); | |
| ae_bool _polynomialfitreport_init_copy(polynomialfitreport* dst, polynomial
fitreport* src, ae_state *_state, ae_bool make_automatic); | | ae_bool _polynomialfitreport_init_copy(polynomialfitreport* dst, polynomial
fitreport* src, ae_state *_state, ae_bool make_automatic); | |
| void _polynomialfitreport_clear(polynomialfitreport* p); | | void _polynomialfitreport_clear(polynomialfitreport* p); | |
|
| | | ae_bool _barycentricfitreport_init(barycentricfitreport* p, ae_state *_stat | |
| | | e, ae_bool make_automatic); | |
| | | ae_bool _barycentricfitreport_init_copy(barycentricfitreport* dst, barycent | |
| | | ricfitreport* src, ae_state *_state, ae_bool make_automatic); | |
| | | void _barycentricfitreport_clear(barycentricfitreport* p); | |
| | | ae_bool _spline1dfitreport_init(spline1dfitreport* p, ae_state *_state, ae_ | |
| | | bool make_automatic); | |
| | | ae_bool _spline1dfitreport_init_copy(spline1dfitreport* dst, spline1dfitrep | |
| | | ort* src, ae_state *_state, ae_bool make_automatic); | |
| | | void _spline1dfitreport_clear(spline1dfitreport* p); | |
| | | ae_bool _lsfitreport_init(lsfitreport* p, ae_state *_state, ae_bool make_au | |
| | | tomatic); | |
| | | ae_bool _lsfitreport_init_copy(lsfitreport* dst, lsfitreport* src, ae_state | |
| | | *_state, ae_bool make_automatic); | |
| | | void _lsfitreport_clear(lsfitreport* p); | |
| | | ae_bool _lsfitstate_init(lsfitstate* p, ae_state *_state, ae_bool make_auto | |
| | | matic); | |
| | | ae_bool _lsfitstate_init_copy(lsfitstate* dst, lsfitstate* src, ae_state *_ | |
| | | state, ae_bool make_automatic); | |
| | | void _lsfitstate_clear(lsfitstate* p); | |
| double barycentriccalc(barycentricinterpolant* b, | | double barycentriccalc(barycentricinterpolant* b, | |
| double t, | | double t, | |
| ae_state *_state); | | ae_state *_state); | |
| void barycentricdiff1(barycentricinterpolant* b, | | void barycentricdiff1(barycentricinterpolant* b, | |
| double t, | | double t, | |
| double* f, | | double* f, | |
| double* df, | | double* df, | |
| ae_state *_state); | | ae_state *_state); | |
| void barycentricdiff2(barycentricinterpolant* b, | | void barycentricdiff2(barycentricinterpolant* b, | |
| double t, | | double t, | |
| | | | |
| skipping to change at line 3227 | | skipping to change at line 3733 | |
|---|
| /* Real */ ae_vector* w, | | /* Real */ ae_vector* w, | |
| ae_int_t n, | | ae_int_t n, | |
| barycentricinterpolant* b, | | barycentricinterpolant* b, | |
| ae_state *_state); | | ae_state *_state); | |
| void barycentricbuildfloaterhormann(/* Real */ ae_vector* x, | | void barycentricbuildfloaterhormann(/* Real */ ae_vector* x, | |
| /* Real */ ae_vector* y, | | /* Real */ ae_vector* y, | |
| ae_int_t n, | | ae_int_t n, | |
| ae_int_t d, | | ae_int_t d, | |
| barycentricinterpolant* b, | | barycentricinterpolant* b, | |
| ae_state *_state); | | ae_state *_state); | |
|
| void barycentricfitfloaterhormannwc(/* Real */ ae_vector* x, | | | |
| /* Real */ ae_vector* y, | | | |
| /* Real */ ae_vector* w, | | | |
| ae_int_t n, | | | |
| /* Real */ ae_vector* xc, | | | |
| /* Real */ ae_vector* yc, | | | |
| /* Integer */ ae_vector* dc, | | | |
| ae_int_t k, | | | |
| ae_int_t m, | | | |
| ae_int_t* info, | | | |
| barycentricinterpolant* b, | | | |
| barycentricfitreport* rep, | | | |
| ae_state *_state); | | | |
| void barycentricfitfloaterhormann(/* Real */ ae_vector* x, | | | |
| /* Real */ ae_vector* y, | | | |
| ae_int_t n, | | | |
| ae_int_t m, | | | |
| ae_int_t* info, | | | |
| barycentricinterpolant* b, | | | |
| barycentricfitreport* rep, | | | |
| ae_state *_state); | | | |
| void barycentriccopy(barycentricinterpolant* b, | | void barycentriccopy(barycentricinterpolant* b, | |
| barycentricinterpolant* b2, | | barycentricinterpolant* b2, | |
| ae_state *_state); | | ae_state *_state); | |
| ae_bool _barycentricinterpolant_init(barycentricinterpolant* p, ae_state *_
state, ae_bool make_automatic); | | ae_bool _barycentricinterpolant_init(barycentricinterpolant* p, ae_state *_
state, ae_bool make_automatic); | |
| ae_bool _barycentricinterpolant_init_copy(barycentricinterpolant* dst, bary
centricinterpolant* src, ae_state *_state, ae_bool make_automatic); | | ae_bool _barycentricinterpolant_init_copy(barycentricinterpolant* dst, bary
centricinterpolant* src, ae_state *_state, ae_bool make_automatic); | |
| void _barycentricinterpolant_clear(barycentricinterpolant* p); | | void _barycentricinterpolant_clear(barycentricinterpolant* p); | |
|
| ae_bool _barycentricfitreport_init(barycentricfitreport* p, ae_state *_stat | | void polynomialbar2cheb(barycentricinterpolant* p, | |
| e, ae_bool make_automatic); | | double a, | |
| ae_bool _barycentricfitreport_init_copy(barycentricfitreport* dst, barycent | | double b, | |
| ricfitreport* src, ae_state *_state, ae_bool make_automatic); | | /* Real */ ae_vector* t, | |
| void _barycentricfitreport_clear(barycentricfitreport* p); | | | |
| void spline2dbuildbilinear(/* Real */ ae_vector* x, | | | |
| /* Real */ ae_vector* y, | | | |
| /* Real */ ae_matrix* f, | | | |
| ae_int_t m, | | | |
| ae_int_t n, | | | |
| spline2dinterpolant* c, | | | |
| ae_state *_state); | | ae_state *_state); | |
|
| void spline2dbuildbicubic(/* Real */ ae_vector* x, | | void polynomialcheb2bar(/* Real */ ae_vector* t, | |
| /* Real */ ae_vector* y, | | | |
| /* Real */ ae_matrix* f, | | | |
| ae_int_t m, | | | |
| ae_int_t n, | | ae_int_t n, | |
|
| spline2dinterpolant* c, | | double a, | |
| | | double b, | |
| | | barycentricinterpolant* p, | |
| ae_state *_state); | | ae_state *_state); | |
|
| double spline2dcalc(spline2dinterpolant* c, | | void polynomialbar2pow(barycentricinterpolant* p, | |
| double x, | | double c, | |
| double y, | | double s, | |
| | | /* Real */ ae_vector* a, | |
| ae_state *_state); | | ae_state *_state); | |
|
| void spline2ddiff(spline2dinterpolant* c, | | void polynomialpow2bar(/* Real */ ae_vector* a, | |
| double x, | | ae_int_t n, | |
| double y, | | double c, | |
| double* f, | | double s, | |
| double* fx, | | barycentricinterpolant* p, | |
| double* fy, | | | |
| double* fxy, | | | |
| ae_state *_state); | | ae_state *_state); | |
|
| void spline2dunpack(spline2dinterpolant* c, | | void polynomialbuild(/* Real */ ae_vector* x, | |
| ae_int_t* m, | | /* Real */ ae_vector* y, | |
| ae_int_t* n, | | ae_int_t n, | |
| /* Real */ ae_matrix* tbl, | | barycentricinterpolant* p, | |
| ae_state *_state); | | ae_state *_state); | |
|
| void spline2dlintransxy(spline2dinterpolant* c, | | void polynomialbuildeqdist(double a, | |
| double ax, | | double b, | |
| double bx, | | /* Real */ ae_vector* y, | |
| double ay, | | ae_int_t n, | |
| double by, | | barycentricinterpolant* p, | |
| ae_state *_state); | | ae_state *_state); | |
|
| void spline2dlintransf(spline2dinterpolant* c, | | void polynomialbuildcheb1(double a, | |
| double a, | | | |
| double b, | | double b, | |
|
| | | /* Real */ ae_vector* y, | |
| | | ae_int_t n, | |
| | | barycentricinterpolant* p, | |
| ae_state *_state); | | ae_state *_state); | |
|
| void spline2dcopy(spline2dinterpolant* c, | | void polynomialbuildcheb2(double a, | |
| spline2dinterpolant* cc, | | double b, | |
| | | /* Real */ ae_vector* y, | |
| | | ae_int_t n, | |
| | | barycentricinterpolant* p, | |
| ae_state *_state); | | ae_state *_state); | |
|
| void spline2dresamplebicubic(/* Real */ ae_matrix* a, | | double polynomialcalceqdist(double a, | |
| ae_int_t oldheight, | | double b, | |
| ae_int_t oldwidth, | | /* Real */ ae_vector* f, | |
| /* Real */ ae_matrix* b, | | ae_int_t n, | |
| ae_int_t newheight, | | double t, | |
| ae_int_t newwidth, | | | |
| ae_state *_state); | | ae_state *_state); | |
|
| void spline2dresamplebilinear(/* Real */ ae_matrix* a, | | double polynomialcalccheb1(double a, | |
| ae_int_t oldheight, | | double b, | |
| ae_int_t oldwidth, | | /* Real */ ae_vector* f, | |
| /* Real */ ae_matrix* b, | | ae_int_t n, | |
| ae_int_t newheight, | | double t, | |
| ae_int_t newwidth, | | ae_state *_state); | |
| | | double polynomialcalccheb2(double a, | |
| | | double b, | |
| | | /* Real */ ae_vector* f, | |
| | | ae_int_t n, | |
| | | double t, | |
| ae_state *_state); | | ae_state *_state); | |
|
| ae_bool _spline2dinterpolant_init(spline2dinterpolant* p, ae_state *_state, | | | |
| ae_bool make_automatic); | | | |
| ae_bool _spline2dinterpolant_init_copy(spline2dinterpolant* dst, spline2din | | | |
| terpolant* src, ae_state *_state, ae_bool make_automatic); | | | |
| void _spline2dinterpolant_clear(spline2dinterpolant* p); | | | |
| void spline1dbuildlinear(/* Real */ ae_vector* x, | | void spline1dbuildlinear(/* Real */ ae_vector* x, | |
| /* Real */ ae_vector* y, | | /* Real */ ae_vector* y, | |
| ae_int_t n, | | ae_int_t n, | |
| spline1dinterpolant* c, | | spline1dinterpolant* c, | |
| ae_state *_state); | | ae_state *_state); | |
| void spline1dbuildcubic(/* Real */ ae_vector* x, | | void spline1dbuildcubic(/* Real */ ae_vector* x, | |
| /* Real */ ae_vector* y, | | /* Real */ ae_vector* y, | |
| ae_int_t n, | | ae_int_t n, | |
| ae_int_t boundltype, | | ae_int_t boundltype, | |
| double boundl, | | double boundl, | |
| | | | |
| skipping to change at line 3405 | | skipping to change at line 3889 | |
|---|
| /* Real */ ae_vector* y, | | /* Real */ ae_vector* y, | |
| /* Real */ ae_vector* d, | | /* Real */ ae_vector* d, | |
| ae_int_t n, | | ae_int_t n, | |
| spline1dinterpolant* c, | | spline1dinterpolant* c, | |
| ae_state *_state); | | ae_state *_state); | |
| void spline1dbuildakima(/* Real */ ae_vector* x, | | void spline1dbuildakima(/* Real */ ae_vector* x, | |
| /* Real */ ae_vector* y, | | /* Real */ ae_vector* y, | |
| ae_int_t n, | | ae_int_t n, | |
| spline1dinterpolant* c, | | spline1dinterpolant* c, | |
| ae_state *_state); | | ae_state *_state); | |
|
| void spline1dfitcubicwc(/* Real */ ae_vector* x, | | | |
| /* Real */ ae_vector* y, | | | |
| /* Real */ ae_vector* w, | | | |
| ae_int_t n, | | | |
| /* Real */ ae_vector* xc, | | | |
| /* Real */ ae_vector* yc, | | | |
| /* Integer */ ae_vector* dc, | | | |
| ae_int_t k, | | | |
| ae_int_t m, | | | |
| ae_int_t* info, | | | |
| spline1dinterpolant* s, | | | |
| spline1dfitreport* rep, | | | |
| ae_state *_state); | | | |
| void spline1dfithermitewc(/* Real */ ae_vector* x, | | | |
| /* Real */ ae_vector* y, | | | |
| /* Real */ ae_vector* w, | | | |
| ae_int_t n, | | | |
| /* Real */ ae_vector* xc, | | | |
| /* Real */ ae_vector* yc, | | | |
| /* Integer */ ae_vector* dc, | | | |
| ae_int_t k, | | | |
| ae_int_t m, | | | |
| ae_int_t* info, | | | |
| spline1dinterpolant* s, | | | |
| spline1dfitreport* rep, | | | |
| ae_state *_state); | | | |
| void spline1dfitcubic(/* Real */ ae_vector* x, | | | |
| /* Real */ ae_vector* y, | | | |
| ae_int_t n, | | | |
| ae_int_t m, | | | |
| ae_int_t* info, | | | |
| spline1dinterpolant* s, | | | |
| spline1dfitreport* rep, | | | |
| ae_state *_state); | | | |
| void spline1dfithermite(/* Real */ ae_vector* x, | | | |
| /* Real */ ae_vector* y, | | | |
| ae_int_t n, | | | |
| ae_int_t m, | | | |
| ae_int_t* info, | | | |
| spline1dinterpolant* s, | | | |
| spline1dfitreport* rep, | | | |
| ae_state *_state); | | | |
| double spline1dcalc(spline1dinterpolant* c, double x, ae_state *_state); | | double spline1dcalc(spline1dinterpolant* c, double x, ae_state *_state); | |
| void spline1ddiff(spline1dinterpolant* c, | | void spline1ddiff(spline1dinterpolant* c, | |
| double x, | | double x, | |
| double* s, | | double* s, | |
| double* ds, | | double* ds, | |
| double* d2s, | | double* d2s, | |
| ae_state *_state); | | ae_state *_state); | |
| void spline1dcopy(spline1dinterpolant* c, | | void spline1dcopy(spline1dinterpolant* c, | |
| spline1dinterpolant* cc, | | spline1dinterpolant* cc, | |
| ae_state *_state); | | ae_state *_state); | |
| | | | |
| skipping to change at line 3472 | | skipping to change at line 3914 | |
|---|
| double a, | | double a, | |
| double b, | | double b, | |
| ae_state *_state); | | ae_state *_state); | |
| void spline1dlintransy(spline1dinterpolant* c, | | void spline1dlintransy(spline1dinterpolant* c, | |
| double a, | | double a, | |
| double b, | | double b, | |
| ae_state *_state); | | ae_state *_state); | |
| double spline1dintegrate(spline1dinterpolant* c, | | double spline1dintegrate(spline1dinterpolant* c, | |
| double x, | | double x, | |
| ae_state *_state); | | ae_state *_state); | |
|
| | | void spline1dconvdiffinternal(/* Real */ ae_vector* xold, | |
| | | /* Real */ ae_vector* yold, | |
| | | /* Real */ ae_vector* dold, | |
| | | ae_int_t n, | |
| | | /* Real */ ae_vector* x2, | |
| | | ae_int_t n2, | |
| | | /* Real */ ae_vector* y, | |
| | | ae_bool needy, | |
| | | /* Real */ ae_vector* d1, | |
| | | ae_bool needd1, | |
| | | /* Real */ ae_vector* d2, | |
| | | ae_bool needd2, | |
| | | ae_state *_state); | |
| | | void heapsortdpoints(/* Real */ ae_vector* x, | |
| | | /* Real */ ae_vector* y, | |
| | | /* Real */ ae_vector* d, | |
| | | ae_int_t n, | |
| | | 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); | |
|
| ae_bool _spline1dfitreport_init(spline1dfitreport* p, ae_state *_state, ae_ | | void spline2dbuildbilinear(/* Real */ ae_vector* x, | |
| bool make_automatic); | | /* Real */ ae_vector* y, | |
| ae_bool _spline1dfitreport_init_copy(spline1dfitreport* dst, spline1dfitrep | | /* Real */ ae_matrix* f, | |
| ort* src, ae_state *_state, ae_bool make_automatic); | | ae_int_t m, | |
| void _spline1dfitreport_clear(spline1dfitreport* p); | | ae_int_t n, | |
| | | spline2dinterpolant* c, | |
| | | ae_state *_state); | |
| | | void spline2dbuildbicubic(/* Real */ ae_vector* x, | |
| | | /* Real */ ae_vector* y, | |
| | | /* Real */ ae_matrix* f, | |
| | | ae_int_t m, | |
| | | ae_int_t n, | |
| | | spline2dinterpolant* c, | |
| | | ae_state *_state); | |
| | | double spline2dcalc(spline2dinterpolant* c, | |
| | | double x, | |
| | | double y, | |
| | | ae_state *_state); | |
| | | void spline2ddiff(spline2dinterpolant* c, | |
| | | double x, | |
| | | double y, | |
| | | double* f, | |
| | | double* fx, | |
| | | double* fy, | |
| | | double* fxy, | |
| | | ae_state *_state); | |
| | | void spline2dunpack(spline2dinterpolant* c, | |
| | | ae_int_t* m, | |
| | | ae_int_t* n, | |
| | | /* Real */ ae_matrix* tbl, | |
| | | ae_state *_state); | |
| | | void spline2dlintransxy(spline2dinterpolant* c, | |
| | | double ax, | |
| | | double bx, | |
| | | double ay, | |
| | | double by, | |
| | | ae_state *_state); | |
| | | void spline2dlintransf(spline2dinterpolant* c, | |
| | | double a, | |
| | | double b, | |
| | | ae_state *_state); | |
| | | void spline2dcopy(spline2dinterpolant* c, | |
| | | spline2dinterpolant* cc, | |
| | | ae_state *_state); | |
| | | void spline2dresamplebicubic(/* Real */ ae_matrix* a, | |
| | | ae_int_t oldheight, | |
| | | ae_int_t oldwidth, | |
| | | /* Real */ ae_matrix* b, | |
| | | ae_int_t newheight, | |
| | | ae_int_t newwidth, | |
| | | ae_state *_state); | |
| | | void spline2dresamplebilinear(/* Real */ ae_matrix* a, | |
| | | ae_int_t oldheight, | |
| | | ae_int_t oldwidth, | |
| | | /* Real */ ae_matrix* b, | |
| | | ae_int_t newheight, | |
| | | ae_int_t newwidth, | |
| | | ae_state *_state); | |
| | | ae_bool _spline2dinterpolant_init(spline2dinterpolant* p, ae_state *_state, | |
| | | ae_bool make_automatic); | |
| | | ae_bool _spline2dinterpolant_init_copy(spline2dinterpolant* dst, spline2din | |
| | | terpolant* src, ae_state *_state, ae_bool make_automatic); | |
| | | void _spline2dinterpolant_clear(spline2dinterpolant* p); | |
| double idwcalc(idwinterpolant* z, | | double idwcalc(idwinterpolant* z, | |
| /* Real */ ae_vector* x, | | /* Real */ ae_vector* x, | |
| ae_state *_state); | | ae_state *_state); | |
| void idwbuildmodifiedshepard(/* Real */ ae_matrix* xy, | | void idwbuildmodifiedshepard(/* Real */ ae_matrix* xy, | |
| ae_int_t n, | | ae_int_t n, | |
| ae_int_t nx, | | ae_int_t nx, | |
| ae_int_t d, | | ae_int_t d, | |
| ae_int_t nq, | | ae_int_t nq, | |
| ae_int_t nw, | | ae_int_t nw, | |
| idwinterpolant* z, | | idwinterpolant* z, | |
| | | | |
| End of changes. 161 change blocks. |
|---|
| 1426 lines changed or deleted | | 1969 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 "solvers.h" | | #include "alglibmisc.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 58 | | skipping to change at line 57 | |
|---|
| ae_int_t q; | | ae_int_t q; | |
| ae_int_t p; | | ae_int_t p; | |
| ae_vector rho; | | ae_vector rho; | |
| ae_matrix y; | | ae_matrix y; | |
| ae_matrix s; | | ae_matrix s; | |
| ae_vector theta; | | ae_vector theta; | |
| ae_vector d; | | ae_vector d; | |
| double stp; | | double stp; | |
| ae_vector work; | | ae_vector work; | |
| double fold; | | double fold; | |
|
| | | ae_int_t prectype; | |
| double gammak; | | double gammak; | |
|
| | | ae_matrix denseh; | |
| | | ae_vector autobuf; | |
| ae_vector x; | | ae_vector x; | |
| double f; | | double f; | |
| ae_vector g; | | ae_vector g; | |
| ae_bool needfg; | | ae_bool needfg; | |
| ae_bool xupdated; | | ae_bool xupdated; | |
| rcommstate rstate; | | rcommstate rstate; | |
| ae_int_t repiterationscount; | | ae_int_t repiterationscount; | |
| ae_int_t repnfev; | | ae_int_t repnfev; | |
| ae_int_t repterminationtype; | | ae_int_t repterminationtype; | |
| linminstate lstate; | | linminstate lstate; | |
| } minlbfgsstate; | | } minlbfgsstate; | |
| typedef struct | | typedef struct | |
| { | | { | |
| 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_bool wrongparams; | | | |
| ae_int_t n; | | ae_int_t n; | |
| ae_int_t m; | | ae_int_t m; | |
|
| | | double diffstep; | |
| 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; | |
| double stpmax; | | double stpmax; | |
|
| ae_int_t flags; | | ae_int_t maxmodelage; | |
| ae_int_t usermode; | | ae_bool makeadditers; | |
| ae_vector x; | | ae_vector x; | |
| double f; | | double f; | |
| ae_vector fi; | | ae_vector fi; | |
| ae_matrix j; | | ae_matrix j; | |
| ae_matrix h; | | ae_matrix h; | |
| ae_vector g; | | ae_vector g; | |
| ae_bool needf; | | ae_bool needf; | |
| ae_bool needfg; | | ae_bool needfg; | |
| ae_bool needfgh; | | ae_bool needfgh; | |
| ae_bool needfij; | | ae_bool needfij; | |
|
| | | ae_bool needfi; | |
| ae_bool xupdated; | | ae_bool xupdated; | |
|
| minlbfgsstate internalstate; | | ae_int_t algomode; | |
| minlbfgsreport internalrep; | | ae_bool hasf; | |
| ae_vector xprec; | | ae_bool hasfi; | |
| | | ae_bool hasg; | |
| ae_vector xbase; | | ae_vector xbase; | |
|
| ae_vector xdir; | | double fbase; | |
| | | ae_vector fibase; | |
| ae_vector gbase; | | ae_vector gbase; | |
|
| ae_vector xprev; | | ae_matrix quadraticmodel; | |
| double fprev; | | double lambdav; | |
| ae_matrix rawmodel; | | double nu; | |
| ae_matrix model; | | ae_matrix dampedmodel; | |
| ae_vector work; | | ae_int_t modelage; | |
| rcommstate rstate; | | ae_vector xdir; | |
| | | ae_vector deltax; | |
| | | ae_vector deltaf; | |
| | | ae_bool deltaxready; | |
| | | ae_bool deltafready; | |
| ae_int_t repiterationscount; | | ae_int_t repiterationscount; | |
| ae_int_t repterminationtype; | | ae_int_t repterminationtype; | |
| ae_int_t repnfunc; | | ae_int_t repnfunc; | |
| ae_int_t repnjac; | | ae_int_t repnjac; | |
| ae_int_t repngrad; | | ae_int_t repngrad; | |
| ae_int_t repnhess; | | ae_int_t repnhess; | |
| ae_int_t repncholesky; | | ae_int_t repncholesky; | |
|
| ae_int_t solverinfo; | | rcommstate rstate; | |
| densesolverreport solverrep; | | ae_vector choleskybuf; | |
| ae_int_t invinfo; | | double actualdecrease; | |
| matinvreport invrep; | | double predicteddecrease; | |
| | | ae_vector fm2; | |
| | | ae_vector fm1; | |
| | | ae_vector fp2; | |
| | | ae_vector fp1; | |
| | | minlbfgsstate internalstate; | |
| | | minlbfgsreport internalrep; | |
| } minlmstate; | | } minlmstate; | |
| typedef struct | | typedef struct | |
| { | | { | |
| ae_int_t iterationscount; | | ae_int_t iterationscount; | |
| ae_int_t terminationtype; | | ae_int_t terminationtype; | |
| ae_int_t nfunc; | | ae_int_t nfunc; | |
| ae_int_t njac; | | ae_int_t njac; | |
| ae_int_t ngrad; | | ae_int_t ngrad; | |
| ae_int_t nhess; | | ae_int_t nhess; | |
| ae_int_t ncholesky; | | ae_int_t ncholesky; | |
| | | | |
| skipping to change at line 301 | | skipping to change at line 316 | |
|---|
| minlbfgsreport(const minlbfgsreport &rhs); | | minlbfgsreport(const minlbfgsreport &rhs); | |
| minlbfgsreport& operator=(const minlbfgsreport &rhs); | | minlbfgsreport& operator=(const minlbfgsreport &rhs); | |
| virtual ~minlbfgsreport(); | | virtual ~minlbfgsreport(); | |
| ae_int_t &iterationscount; | | ae_int_t &iterationscount; | |
| ae_int_t &nfev; | | ae_int_t &nfev; | |
| ae_int_t &terminationtype; | | ae_int_t &terminationtype; | |
| | | | |
| }; | | }; | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
|
| | | Levenberg-Marquardt optimizer. | |
| | | | |
|
| | | This structure should be created using one of the MinLMCreate???() | |
| | | functions. You should not access its fields directly; use ALGLIB functions | |
| | | to work with it. | |
| *************************************************************************/ | | *************************************************************************/ | |
| class _minlmstate_owner | | class _minlmstate_owner | |
| { | | { | |
| public: | | public: | |
| _minlmstate_owner(); | | _minlmstate_owner(); | |
| _minlmstate_owner(const _minlmstate_owner &rhs); | | _minlmstate_owner(const _minlmstate_owner &rhs); | |
| _minlmstate_owner& operator=(const _minlmstate_owner &rhs); | | _minlmstate_owner& operator=(const _minlmstate_owner &rhs); | |
| virtual ~_minlmstate_owner(); | | virtual ~_minlmstate_owner(); | |
| alglib_impl::minlmstate* c_ptr(); | | alglib_impl::minlmstate* c_ptr(); | |
| alglib_impl::minlmstate* c_ptr() const; | | alglib_impl::minlmstate* c_ptr() const; | |
| | | | |
| skipping to change at line 325 | | skipping to change at line 344 | |
|---|
| class minlmstate : public _minlmstate_owner | | class minlmstate : public _minlmstate_owner | |
| { | | { | |
| public: | | public: | |
| minlmstate(); | | minlmstate(); | |
| minlmstate(const minlmstate &rhs); | | minlmstate(const minlmstate &rhs); | |
| minlmstate& operator=(const minlmstate &rhs); | | minlmstate& operator=(const minlmstate &rhs); | |
| virtual ~minlmstate(); | | virtual ~minlmstate(); | |
| ae_bool &needf; | | ae_bool &needf; | |
| ae_bool &needfg; | | ae_bool &needfg; | |
| ae_bool &needfgh; | | ae_bool &needfgh; | |
|
| | | ae_bool &needfi; | |
| ae_bool &needfij; | | ae_bool &needfij; | |
| ae_bool &xupdated; | | ae_bool &xupdated; | |
| double &f; | | double &f; | |
| real_1d_array fi; | | real_1d_array fi; | |
| real_1d_array g; | | real_1d_array g; | |
| real_2d_array h; | | real_2d_array h; | |
| real_2d_array j; | | real_2d_array j; | |
| real_1d_array x; | | real_1d_array x; | |
| | | | |
| }; | | }; | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
|
| | | Optimization report, filled by MinLMResults() function | |
| | | | |
|
| | | FIELDS: | |
| | | * TerminationType, completetion code: | |
| | | * -9 derivative correctness check failed; | |
| | | see Rep.WrongNum, Rep.WrongI, Rep.WrongJ for | |
| | | more information. | |
| | | * 1 relative function improvement is no more than | |
| | | EpsF. | |
| | | * 2 relative step is no more than EpsX. | |
| | | * 4 gradient is no more than EpsG. | |
| | | * 5 MaxIts steps was taken | |
| | | * 7 stopping conditions are too stringent, | |
| | | further improvement is impossible | |
| | | * IterationsCount, contains iterations count | |
| | | * NFunc, number of function calculations | |
| | | * NJac, number of Jacobi matrix calculations | |
| | | * NGrad, number of gradient calculations | |
| | | * NHess, number of Hessian calculations | |
| | | * NCholesky, number of Cholesky decomposition calculations | |
| *************************************************************************/ | | *************************************************************************/ | |
| class _minlmreport_owner | | class _minlmreport_owner | |
| { | | { | |
| public: | | public: | |
| _minlmreport_owner(); | | _minlmreport_owner(); | |
| _minlmreport_owner(const _minlmreport_owner &rhs); | | _minlmreport_owner(const _minlmreport_owner &rhs); | |
| _minlmreport_owner& operator=(const _minlmreport_owner &rhs); | | _minlmreport_owner& operator=(const _minlmreport_owner &rhs); | |
| virtual ~_minlmreport_owner(); | | virtual ~_minlmreport_owner(); | |
| alglib_impl::minlmreport* c_ptr(); | | alglib_impl::minlmreport* c_ptr(); | |
| alglib_impl::minlmreport* c_ptr() const; | | alglib_impl::minlmreport* c_ptr() const; | |
| | | | |
| skipping to change at line 599 | | skipping to change at line 638 | |
|---|
| large steps which leads to overflow. This function allows us to reject | | large steps which leads to overflow. This function allows us to reject | |
| steps that are too large (and therefore expose us to the possible | | steps that are too large (and therefore expose us to the possible | |
| overflow) without actually calculating function value at the x+stp*d. | | overflow) without actually calculating function value at the x+stp*d. | |
| | | | |
| -- ALGLIB -- | | -- ALGLIB -- | |
| Copyright 02.04.2010 by Bochkanov Sergey | | Copyright 02.04.2010 by Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
| void minlbfgssetstpmax(const minlbfgsstate &state, const double stpmax); | | void minlbfgssetstpmax(const minlbfgsstate &state, const double stpmax); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
|
| | | Modification of the preconditioner: | |
| | | default preconditioner (simple scaling) is used. | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | State - structure which stores algorithm state | |
| | | | |
| | | After call to this function preconditioner is changed to the default one. | |
| | | | |
| | | NOTE: you can change preconditioner "on the fly", during algorithm | |
| | | iterations. | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 13.10.2010 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void minlbfgssetdefaultpreconditioner(const minlbfgsstate &state); | |
| | | | |
| | | /************************************************************************* | |
| | | Modification of the preconditioner: | |
| | | Cholesky factorization of approximate Hessian is used. | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | State - structure which stores algorithm state | |
| | | P - triangular preconditioner, Cholesky factorization of | |
| | | the approximate Hessian. array[0..N-1,0..N-1], | |
| | | (if larger, only leading N elements are used). | |
| | | IsUpper - whether upper or lower triangle of P is given | |
| | | (other triangle is not referenced) | |
| | | | |
| | | After call to this function preconditioner is changed to P (P is copied | |
| | | into the internal buffer). | |
| | | | |
| | | NOTE: you can change preconditioner "on the fly", during algorithm | |
| | | iterations. | |
| | | | |
| | | NOTE 2: P should be nonsingular. Exception will be thrown otherwise. It | |
| | | also should be well conditioned, although only strict non-singularity is | |
| | | tested. | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 13.10.2010 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void minlbfgssetcholeskypreconditioner(const minlbfgsstate &state, const re | |
| | | al_2d_array &p, const bool isupper); | |
| | | | |
| | | /************************************************************************* | |
| This function provides reverse communication interface | | This function provides reverse communication interface | |
| Reverse communication interface is not documented or recommended to use. | | Reverse communication interface is not documented or recommended to use. | |
| See below for functions which provide better documented API | | See below for functions which provide better documented API | |
| *************************************************************************/ | | *************************************************************************/ | |
| bool minlbfgsiteration(const minlbfgsstate &state); | | bool minlbfgsiteration(const minlbfgsstate &state); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
| This family of functions is used to launcn iterations of nonlinear optimize
r | | This family of functions is used to launcn iterations of nonlinear optimize
r | |
| | | | |
| These functions accept following parameters: | | These functions accept following parameters: | |
| | | | |
| skipping to change at line 684 | | skipping to change at line 767 | |
|---|
| INPUT PARAMETERS: | | INPUT PARAMETERS: | |
| State - structure used to store algorithm state | | State - structure used to store algorithm state | |
| X - new starting point. | | X - new starting point. | |
| | | | |
| -- ALGLIB -- | | -- ALGLIB -- | |
| Copyright 30.07.2010 by Bochkanov Sergey | | Copyright 30.07.2010 by Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
| void minlbfgsrestartfrom(const minlbfgsstate &state, const real_1d_array &x
); | | void minlbfgsrestartfrom(const minlbfgsstate &state, const real_1d_array &x
); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
|
| LEVENBERG-MARQUARDT-LIKE METHOD FOR NON-LINEAR OPTIMIZATION | | IMPROVED LEVENBERG-MARQUARDT METHOD FOR | |
| | | NON-LINEAR LEAST SQUARES OPTIMIZATION | |
| | | | |
| DESCRIPTION: | | DESCRIPTION: | |
|
| This function is used to find minimum of general form (not "sum-of- | | This function is used to find minimum of function which is represented as | |
| -squares") function | | sum of squares: | |
| F = F(x[0], ..., x[n-1]) | | F(x) = f[0]^2(x[0],...,x[n-1]) + ... + f[m-1]^2(x[0],...,x[n-1]) | |
| using its gradient and Hessian. Levenberg-Marquardt modification with | | using value of function vector f[] and Jacobian of f[]. | |
| L-BFGS pre-optimization and internal pre-conditioned L-BFGS optimization | | | |
| after each Levenberg-Marquardt step is used. | | | |
| | | | |
| REQUIREMENTS: | | REQUIREMENTS: | |
| This algorithm will request following information during its operation: | | This algorithm will request following information during its operation: | |
| | | | |
|
| * function value F at given point X | | * function vector f[] at given point X | |
| * F and gradient G (simultaneously) at given point X | | * function vector f[] and Jacobian of f[] (simultaneously) at given point | |
| * F, G and Hessian H (simultaneously) at given point X | | | |
| | | | |
| There are several overloaded versions of MinLMOptimize() function which | | There are several overloaded versions of MinLMOptimize() function which | |
| correspond to different LM-like optimization algorithms provided by this | | correspond to different LM-like optimization algorithms provided by this | |
|
| unit. You should choose version which accepts func(), grad() and hess() | | unit. You should choose version which accepts fvec() and jac() callbacks. | |
| function pointers. First pointer is used to calculate F at given point, | | First one is used to calculate f[] at given point, second one calculates | |
| second one calculates F(x) and grad F(x), third one calculates F(x), | | f[] and Jacobian df[i]/dx[j]. | |
| grad F(x), hess F(x). | | | |
| | | | |
|
| You can try to initialize MinLMState structure with FGH-function and then | | You can try to initialize MinLMState structure with VJ function and then | |
| use incorrect version of MinLMOptimize() (for example, version which does | | use incorrect version of MinLMOptimize() (for example, version which | |
| not provide Hessian matrix), but it will lead to exception being thrown | | works with general form function and does not provide Jacobian), but it | |
| after first attempt to calculate Hessian. | | will lead to exception being thrown after first attempt to calculate | |
| | | Jacobian. | |
| | | | |
| USAGE: | | USAGE: | |
|
| 1. User initializes algorithm state with MinLMCreateFGH() call | | 1. User initializes algorithm state with MinLMCreateVJ() call | |
| 2. User tunes solver parameters with MinLMSetCond(), MinLMSetStpMax() and | | 2. User tunes solver parameters with MinLMSetCond(), MinLMSetStpMax() and | |
| other functions | | other functions | |
| 3. User calls MinLMOptimize() function which takes algorithm state and | | 3. User calls MinLMOptimize() function which takes algorithm state and | |
|
| pointers (delegates, etc.) to callback functions. | | callback functions. | |
| 4. User calls MinLMResults() to get solution | | 4. User calls MinLMResults() to get solution | |
| 5. Optionally, user may call MinLMRestartFrom() to solve another problem | | 5. Optionally, user may call MinLMRestartFrom() to solve another problem | |
|
| with same N but another starting point and/or another function. | | with same N/M but another starting point and/or another function. | |
| MinLMRestartFrom() allows to reuse already initialized structure. | | MinLMRestartFrom() allows to reuse already initialized structure. | |
| | | | |
| INPUT PARAMETERS: | | INPUT PARAMETERS: | |
| N - dimension, N>1 | | N - dimension, N>1 | |
| * if given, only leading N elements of X are used | | * if given, only leading N elements of X are used | |
| * if not given, automatically determined from size of X | | * if not given, automatically determined from size of X | |
|
| | | M - number of functions f[i] | |
| X - initial solution, array[0..N-1] | | X - initial solution, array[0..N-1] | |
| | | | |
| OUTPUT PARAMETERS: | | OUTPUT PARAMETERS: | |
| State - structure which stores algorithm state | | State - structure which stores algorithm state | |
| | | | |
| NOTES: | | NOTES: | |
| 1. you may tune stopping conditions with MinLMSetCond() function | | 1. you may tune stopping conditions with MinLMSetCond() function | |
| 2. if target function contains exp() or other fast growing functions, and | | 2. if target function contains exp() or other fast growing functions, and | |
| optimization algorithm makes too large steps which leads to overflow, | | optimization algorithm makes too large steps which leads to overflow, | |
| use MinLMSetStpMax() function to bound algorithm's steps. | | use MinLMSetStpMax() function to bound algorithm's steps. | |
| | | | |
| -- ALGLIB -- | | -- ALGLIB -- | |
| Copyright 30.03.2009 by Bochkanov Sergey | | Copyright 30.03.2009 by Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
|
| void minlmcreatefgh(const ae_int_t n, const real_1d_array &x, minlmstate &s | | void minlmcreatevj(const ae_int_t n, const ae_int_t m, const real_1d_array | |
| tate); | | &x, minlmstate &state); | |
| void minlmcreatefgh(const real_1d_array &x, minlmstate &state); | | void minlmcreatevj(const ae_int_t m, const real_1d_array &x, minlmstate &st | |
| | | ate); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
|
| LEVENBERG-MARQUARDT-LIKE METHOD FOR NON-LINEAR OPTIMIZATION | | IMPROVED LEVENBERG-MARQUARDT METHOD FOR | |
| | | NON-LINEAR LEAST SQUARES OPTIMIZATION | |
| | | | |
| DESCRIPTION: | | DESCRIPTION: | |
| This function is used to find minimum of function which is represented as | | This function is used to find minimum of function which is represented as | |
| sum of squares: | | sum of squares: | |
| F(x) = f[0]^2(x[0],...,x[n-1]) + ... + f[m-1]^2(x[0],...,x[n-1]) | | F(x) = f[0]^2(x[0],...,x[n-1]) + ... + f[m-1]^2(x[0],...,x[n-1]) | |
|
| using value of F(), gradient of F() function vector f[] and Jacobian of | | using value of function vector f[] only. Finite differences are used to | |
| f[]. Levenberg-Marquardt modification with L-BFGS pre-optimization and | | calculate Jacobian. | |
| internal preconditioned L-BFGS optimization after each Levenberg-Marquardt | | | |
| step is used. | | REQUIREMENTS: | |
| | | This algorithm will request following information during its operation: | |
| | | * function vector f[] at given point X | |
| | | | |
| | | There are several overloaded versions of MinLMOptimize() function which | |
| | | correspond to different LM-like optimization algorithms provided by this | |
| | | unit. You should choose version which accepts fvec() callback. | |
| | | | |
| | | You can try to initialize MinLMState structure with VJ function and then | |
| | | use incorrect version of MinLMOptimize() (for example, version which | |
| | | works with general form function and does not accept function vector), but | |
| | | it will lead to exception being thrown after first attempt to calculate | |
| | | Jacobian. | |
| | | | |
| | | USAGE: | |
| | | 1. User initializes algorithm state with MinLMCreateV() call | |
| | | 2. User tunes solver parameters with MinLMSetCond(), MinLMSetStpMax() and | |
| | | other functions | |
| | | 3. User calls MinLMOptimize() function which takes algorithm state and | |
| | | callback functions. | |
| | | 4. User calls MinLMResults() to get solution | |
| | | 5. Optionally, user may call MinLMRestartFrom() to solve another problem | |
| | | with same N/M but another starting point and/or another function. | |
| | | MinLMRestartFrom() allows to reuse already initialized structure. | |
| | | | |
| | | INPUT PARAMETERS: | |
| | | N - dimension, N>1 | |
| | | * if given, only leading N elements of X are used | |
| | | * if not given, automatically determined from size of X | |
| | | M - number of functions f[i] | |
| | | X - initial solution, array[0..N-1] | |
| | | DiffStep- differentiation step, >0 | |
| | | | |
| | | OUTPUT PARAMETERS: | |
| | | State - structure which stores algorithm state | |
| | | | |
| | | See also MinLMIteration, MinLMResults. | |
| | | | |
| | | NOTES: | |
| | | 1. you may tune stopping conditions with MinLMSetCond() function | |
| | | 2. if target function contains exp() or other fast growing functions, and | |
| | | optimization algorithm makes too large steps which leads to overflow, | |
| | | use MinLMSetStpMax() function to bound algorithm's steps. | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 30.03.2009 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void minlmcreatev(const ae_int_t n, const ae_int_t m, const real_1d_array & | |
| | | x, const double diffstep, minlmstate &state); | |
| | | void minlmcreatev(const ae_int_t m, const real_1d_array &x, const double di | |
| | | ffstep, minlmstate &state); | |
| | | | |
| | | /************************************************************************* | |
| | | LEVENBERG-MARQUARDT-LIKE METHOD FOR NON-LINEAR OPTIMIZATION | |
| | | | |
| | | DESCRIPTION: | |
| | | This function is used to find minimum of general form (not "sum-of- | |
| | | -squares") function | |
| | | F = F(x[0], ..., x[n-1]) | |
| | | using its gradient and Hessian. Levenberg-Marquardt modification with | |
| | | L-BFGS pre-optimization and internal pre-conditioned L-BFGS optimization | |
| | | after each Levenberg-Marquardt step is used. | |
| | | | |
| REQUIREMENTS: | | REQUIREMENTS: | |
| This algorithm will request following information during its operation: | | This algorithm will request following information during its operation: | |
| | | | |
| * function value F at given point X | | * function value F at given point X | |
| * F and gradient G (simultaneously) at given point X | | * F and gradient G (simultaneously) at given point X | |
|
| * function vector f[] and Jacobian of f[] (simultaneously) at given point | | * F, G and Hessian H (simultaneously) at given point X | |
| | | | |
| There are several overloaded versions of MinLMOptimize() function which | | There are several overloaded versions of MinLMOptimize() function which | |
| correspond to different LM-like optimization algorithms provided by this | | correspond to different LM-like optimization algorithms provided by this | |
|
| unit. You should choose version which accepts func(), grad() and jac() | | unit. You should choose version which accepts func(), grad() and hess() | |
| function pointers. First pointer is used to calculate F at given point, | | function pointers. First pointer is used to calculate F at given point, | |
|
| second one calculates F(x) and grad F(x), third one calculates f[] and | | second one calculates F(x) and grad F(x), third one calculates F(x), | |
| Jacobian df[i]/dx[j]. | | grad F(x), hess F(x). | |
| | | | |
|
| You can try to initialize MinLMState structure with FGJ function and then | | You can try to initialize MinLMState structure with FGH-function and then | |
| use incorrect version of MinLMOptimize() (for example, version which does | | use incorrect version of MinLMOptimize() (for example, version which does | |
|
| not provide gradient), but it will lead to exception being thrown after | | not provide Hessian matrix), but it will lead to exception being thrown | |
| first attempt to calculate gradient. | | after first attempt to calculate Hessian. | |
| | | | |
| USAGE: | | USAGE: | |
|
| 1. User initializes algorithm state with MinLMCreateFGJ() call | | 1. User initializes algorithm state with MinLMCreateFGH() call | |
| 2. User tunes solver parameters with MinLMSetCond(), MinLMSetStpMax() and | | 2. User tunes solver parameters with MinLMSetCond(), MinLMSetStpMax() and | |
| other functions | | other functions | |
| 3. User calls MinLMOptimize() function which takes algorithm state and | | 3. User calls MinLMOptimize() function which takes algorithm state and | |
| pointers (delegates, etc.) to callback functions. | | pointers (delegates, etc.) to callback functions. | |
| 4. User calls MinLMResults() to get solution | | 4. User calls MinLMResults() to get solution | |
| 5. Optionally, user may call MinLMRestartFrom() to solve another problem | | 5. Optionally, user may call MinLMRestartFrom() to solve another problem | |
|
| with same N/M but another starting point and/or another function. | | with same N but another starting point and/or another function. | |
| MinLMRestartFrom() allows to reuse already initialized structure. | | MinLMRestartFrom() allows to reuse already initialized structure. | |
| | | | |
| INPUT PARAMETERS: | | INPUT PARAMETERS: | |
|
| N - dimension, N>1: | | N - dimension, N>1 | |
| * if given, only leading N elements of X are used | | * if given, only leading N elements of X are used | |
| * if not given, automatically determined from size of X | | * if not given, automatically determined from size of X | |
|
| M - number of functions f[i] | | | |
| X - initial solution, array[0..N-1] | | X - initial solution, array[0..N-1] | |
| | | | |
| OUTPUT PARAMETERS: | | OUTPUT PARAMETERS: | |
| State - structure which stores algorithm state | | State - structure which stores algorithm state | |
| | | | |
|
| See also MinLMIteration, MinLMResults. | | | |
| | | | |
| NOTES: | | NOTES: | |
| 1. you may tune stopping conditions with MinLMSetCond() function | | 1. you may tune stopping conditions with MinLMSetCond() function | |
| 2. if target function contains exp() or other fast growing functions, and | | 2. if target function contains exp() or other fast growing functions, and | |
| optimization algorithm makes too large steps which leads to overflow, | | optimization algorithm makes too large steps which leads to overflow, | |
| use MinLMSetStpMax() function to bound algorithm's steps. | | use MinLMSetStpMax() function to bound algorithm's steps. | |
| | | | |
| -- ALGLIB -- | | -- ALGLIB -- | |
| Copyright 30.03.2009 by Bochkanov Sergey | | Copyright 30.03.2009 by Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
|
| void minlmcreatefgj(const ae_int_t n, const ae_int_t m, const real_1d_array | | void minlmcreatefgh(const ae_int_t n, const real_1d_array &x, minlmstate &s | |
| &x, minlmstate &state); | | tate); | |
| void minlmcreatefgj(const ae_int_t m, const real_1d_array &x, minlmstate &s | | void minlmcreatefgh(const real_1d_array &x, minlmstate &state); | |
| tate); | | | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
|
| CLASSIC LEVENBERG-MARQUARDT METHOD FOR NON-LINEAR OPTIMIZATION | | IMPROVED LEVENBERG-MARQUARDT METHOD FOR | |
| | | NON-LINEAR LEAST SQUARES OPTIMIZATION | |
| | | | |
| DESCRIPTION: | | DESCRIPTION: | |
| This function is used to find minimum of function which is represented as | | This function is used to find minimum of function which is represented as | |
| sum of squares: | | sum of squares: | |
| F(x) = f[0]^2(x[0],...,x[n-1]) + ... + f[m-1]^2(x[0],...,x[n-1]) | | F(x) = f[0]^2(x[0],...,x[n-1]) + ... + f[m-1]^2(x[0],...,x[n-1]) | |
|
| using value of F(), function vector f[] and Jacobian of f[]. Classic | | using: | |
| Levenberg-Marquardt method is used. | | * value of function vector f[] | |
| | | * value of Jacobian of f[] | |
| | | * gradient of merit function F(x) | |
| | | | |
| | | This function creates optimizer which uses acceleration strategy 2. Cheap | |
| | | gradient of merit function (which is twice the product of function vector | |
| | | and Jacobian) is used for accelerated iterations (see User Guide for more | |
| | | info on this subject). | |
| | | | |
| REQUIREMENTS: | | REQUIREMENTS: | |
| This algorithm will request following information during its operation: | | This algorithm will request following information during its operation: | |
| | | | |
|
| * function value F at given point X | | * function vector f[] at given point X | |
| * function vector f[] and Jacobian of f[] (simultaneously) at given point | | * function vector f[] and Jacobian of f[] (simultaneously) at given point | |
|
| | | * gradient of | |
| | | | |
| There are several overloaded versions of MinLMOptimize() function which | | There are several overloaded versions of MinLMOptimize() function which | |
| correspond to different LM-like optimization algorithms provided by this | | correspond to different LM-like optimization algorithms provided by this | |
|
| unit. You should choose version which accepts func() and jac() function | | unit. You should choose version which accepts fvec(), jac() and grad() | |
| pointers. First pointer is used to calculate F at given point, second one | | callbacks. First one is used to calculate f[] at given point, second one | |
| calculates calculates f[] and Jacobian df[i]/dx[j]. | | calculates f[] and Jacobian df[i]/dx[j], last one calculates gradient of | |
| | | merit function F(x). | |
| | | | |
|
| You can try to initialize MinLMState structure with FJ function and then | | You can try to initialize MinLMState structure with VJ function and then | |
| use incorrect version of MinLMOptimize() (for example, version which | | use incorrect version of MinLMOptimize() (for example, version which | |
| works with general form function and does not provide Jacobian), but it | | works with general form function and does not provide Jacobian), but it | |
| will lead to exception being thrown after first attempt to calculate | | will lead to exception being thrown after first attempt to calculate | |
| Jacobian. | | Jacobian. | |
| | | | |
| USAGE: | | USAGE: | |
|
| 1. User initializes algorithm state with MinLMCreateFJ() call | | 1. User initializes algorithm state with MinLMCreateVGJ() call | |
| 2. User tunes solver parameters with MinLMSetCond(), MinLMSetStpMax() and | | 2. User tunes solver parameters with MinLMSetCond(), MinLMSetStpMax() and | |
| other functions | | other functions | |
| 3. User calls MinLMOptimize() function which takes algorithm state and | | 3. User calls MinLMOptimize() function which takes algorithm state and | |
|
| pointers (delegates, etc.) to callback functions. | | callback functions. | |
| 4. User calls MinLMResults() to get solution | | 4. User calls MinLMResults() to get solution | |
| 5. Optionally, user may call MinLMRestartFrom() to solve another problem | | 5. Optionally, user may call MinLMRestartFrom() to solve another problem | |
| with same N/M but another starting point and/or another function. | | with same N/M but another starting point and/or another function. | |
| MinLMRestartFrom() allows to reuse already initialized structure. | | MinLMRestartFrom() allows to reuse already initialized structure. | |
| | | | |
| INPUT PARAMETERS: | | INPUT PARAMETERS: | |
| N - dimension, N>1 | | N - dimension, N>1 | |
| * if given, only leading N elements of X are used | | * if given, only leading N elements of X are used | |
| * if not given, automatically determined from size of X | | * if not given, automatically determined from size of X | |
| M - number of functions f[i] | | M - number of functions f[i] | |
| X - initial solution, array[0..N-1] | | X - initial solution, array[0..N-1] | |
| | | | |
| OUTPUT PARAMETERS: | | OUTPUT PARAMETERS: | |
| State - structure which stores algorithm state | | State - structure which stores algorithm state | |
| | | | |
|
| See also MinLMIteration, MinLMResults. | | | |
| | | | |
| NOTES: | | NOTES: | |
| 1. you may tune stopping conditions with MinLMSetCond() function | | 1. you may tune stopping conditions with MinLMSetCond() function | |
| 2. if target function contains exp() or other fast growing functions, and | | 2. if target function contains exp() or other fast growing functions, and | |
| optimization algorithm makes too large steps which leads to overflow, | | optimization algorithm makes too large steps which leads to overflow, | |
| use MinLMSetStpMax() function to bound algorithm's steps. | | use MinLMSetStpMax() function to bound algorithm's steps. | |
| | | | |
| -- ALGLIB -- | | -- ALGLIB -- | |
| Copyright 30.03.2009 by Bochkanov Sergey | | Copyright 30.03.2009 by Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
|
| | | void minlmcreatevgj(const ae_int_t n, const ae_int_t m, const real_1d_array | |
| | | &x, minlmstate &state); | |
| | | void minlmcreatevgj(const ae_int_t m, const real_1d_array &x, minlmstate &s | |
| | | tate); | |
| | | | |
| | | /************************************************************************* | |
| | | LEVENBERG-MARQUARDT-LIKE METHOD FOR | |
| | | NON-LINEAR LEAST SQUARES OPTIMIZATION | |
| | | | |
| | | DESCRIPTION: | |
| | | | |
| | | This function is used to find minimum of function which is represented as | |
| | | sum of squares: | |
| | | F(x) = f[0]^2(x[0],...,x[n-1]) + ... + f[m-1]^2(x[0],...,x[n-1]) | |
| | | using value of F(), gradient of F(), function vector f[] and Jacobian of | |
| | | f[]. | |
| | | | |
| | | This function is considered obsolete since ALGLIB 3.1.0 and is present for | |
| | | backward compatibility only. We recommend to use MinLMCreateVGJ, which | |
| | | provides similar, but more consistent interface. | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 30.03.2009 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void minlmcreatefgj(const ae_int_t n, const ae_int_t m, const real_1d_array | |
| | | &x, minlmstate &state); | |
| | | void minlmcreatefgj(const ae_int_t m, const real_1d_array &x, minlmstate &s | |
| | | tate); | |
| | | | |
| | | /************************************************************************* | |
| | | CLASSIC LEVENBERG-MARQUARDT METHOD FOR NON-LINEAR OPTIMIZATION | |
| | | | |
| | | DESCRIPTION: | |
| | | This function is used to find minimum of function which is represented as | |
| | | sum of squares: | |
| | | F(x) = f[0]^2(x[0],...,x[n-1]) + ... + f[m-1]^2(x[0],...,x[n-1]) | |
| | | using value of F(), function vector f[] and Jacobian of f[]. Classic | |
| | | Levenberg-Marquardt method is used. | |
| | | | |
| | | This function is considered obsolete since ALGLIB 3.1.0 and is present for | |
| | | backward compatibility only. We recommend to use MinLMCreateVJ, which | |
| | | provides similar, but more consistent and feature-rich interface. | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 30.03.2009 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| void minlmcreatefj(const ae_int_t n, const ae_int_t m, const real_1d_array
&x, minlmstate &state); | | void minlmcreatefj(const ae_int_t n, const ae_int_t m, const real_1d_array
&x, minlmstate &state); | |
| void minlmcreatefj(const ae_int_t m, const real_1d_array &x, minlmstate &st
ate); | | void minlmcreatefj(const ae_int_t m, const real_1d_array &x, minlmstate &st
ate); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
| This function sets stopping conditions for Levenberg-Marquardt optimization | | This function sets stopping conditions for Levenberg-Marquardt optimization | |
| algorithm. | | algorithm. | |
| | | | |
| INPUT PARAMETERS: | | INPUT PARAMETERS: | |
| State - structure which stores algorithm state | | State - structure which stores algorithm state | |
| EpsG - >=0 | | EpsG - >=0 | |
| | | | |
| skipping to change at line 945 | | skipping to change at line 1134 | |
|---|
| NOTE: non-zero StpMax leads to moderate performance degradation because | | NOTE: non-zero StpMax leads to moderate performance degradation because | |
| intermediate step of preconditioned L-BFGS optimization is incompatible | | intermediate step of preconditioned L-BFGS optimization is incompatible | |
| with limits on step size. | | with limits on step size. | |
| | | | |
| -- ALGLIB -- | | -- ALGLIB -- | |
| Copyright 02.04.2010 by Bochkanov Sergey | | Copyright 02.04.2010 by Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
| void minlmsetstpmax(const minlmstate &state, const double stpmax); | | void minlmsetstpmax(const minlmstate &state, const double stpmax); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
|
| | | This function is used to change acceleration settings | |
| | | | |
| | | You can choose between three acceleration strategies: | |
| | | * AccType=0, no acceleration. | |
| | | * AccType=1, secant updates are used to update quadratic model after each | |
| | | iteration. After fixed number of iterations (or after model breakdown) | |
| | | we recalculate quadratic model using analytic Jacobian or finite | |
| | | differences. Number of secant-based iterations depends on optimization | |
| | | settings: about 3 iterations - when we have analytic Jacobian, up to 2*N | |
| | | iterations - when we use finite differences to calculate Jacobian. | |
| | | * AccType=2, after quadratic model is built and LM step is made, we use it | |
| | | as preconditioner for several (5-10) iterations of L-BFGS algorithm. | |
| | | | |
| | | AccType=1 is recommended when Jacobian calculation cost is prohibitive | |
| | | high (several Mx1 function vector calculations followed by several NxN | |
| | | Cholesky factorizations are faster than calculation of one M*N Jacobian). | |
| | | It should also be used when we have no Jacobian, because finite difference | |
| | | approximation takes too much time to compute. | |
| | | | |
| | | AccType=2 is recommended when Jacobian is cheap - much more cheaper than | |
| | | one Cholesky factorization. We can reduce number of Cholesky | |
| | | factorizations at the cost of increased number of Jacobian calculations. | |
| | | Sometimes it helps. | |
| | | | |
| | | Table below list optimization protocols (XYZ protocol corresponds to | |
| | | MinLMCreateXYZ) and acceleration types they support (and use by default). | |
| | | | |
| | | ACCELERATION TYPES SUPPORTED BY OPTIMIZATION PROTOCOLS: | |
| | | | |
| | | protocol 0 1 2 comment | |
| | | V + + | |
| | | VJ + + + | |
| | | FGH + + | |
| | | VGJ + + + special protocol, not for widespread use | |
| | | FJ + + obsolete protocol, not recommended | |
| | | FGJ + + obsolete protocol, not recommended | |
| | | | |
| | | DAFAULT VALUES: | |
| | | | |
| | | protocol 0 1 2 comment | |
| | | V x without acceleration it is so slooooooooow | |
| | | VJ x | |
| | | FGH x | |
| | | VGJ x we've implicitly turned (2) by passing gradient | |
| | | FJ x obsolete protocol, not recommended | |
| | | FGJ x obsolete protocol, not recommended | |
| | | | |
| | | NOTE: this function should be called before optimization. Attempt to call | |
| | | it during algorithm iterations may result in unexpected behavior. | |
| | | | |
| | | NOTE: attempt to call this function with unsupported protocol/acceleration | |
| | | combination will result in exception being thrown. | |
| | | | |
| | | -- ALGLIB -- | |
| | | Copyright 14.10.2010 by Bochkanov Sergey | |
| | | *************************************************************************/ | |
| | | void minlmsetacctype(const minlmstate &state, const ae_int_t acctype); | |
| | | | |
| | | /************************************************************************* | |
| This function provides reverse communication interface | | This function provides reverse communication interface | |
| Reverse communication interface is not documented or recommended to use. | | Reverse communication interface is not documented or recommended to use. | |
| See below for functions which provide better documented API | | See below for functions which provide better documented API | |
| *************************************************************************/ | | *************************************************************************/ | |
| bool minlmiteration(const minlmstate &state); | | bool minlmiteration(const minlmstate &state); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
| This family of functions is used to launcn iterations of nonlinear optimize
r | | This family of functions is used to launcn iterations of nonlinear optimize
r | |
| | | | |
| These functions accept following parameters: | | These functions accept following parameters: | |
| state - algorithm state | | state - algorithm state | |
| func - callback which calculates function (or merit function) | | func - callback which calculates function (or merit function) | |
| value func at given point x | | value func at given point x | |
| grad - callback which calculates function (or merit function) | | grad - callback which calculates function (or merit function) | |
| value func and gradient grad at given point x | | value func and gradient grad at given point x | |
| hess - callback which calculates function (or merit function) | | hess - callback which calculates function (or merit function) | |
| value func, gradient grad and Hessian hess at given point x | | value func, gradient grad and Hessian hess at given point x | |
|
| | | fvec - callback which calculates function vector fi[] | |
| | | at given point x | |
| jac - callback which calculates function vector fi[] | | jac - callback which calculates function vector fi[] | |
| and Jacobian jac at given point x | | and Jacobian jac at given point x | |
| rep - optional callback which is called after each iteration | | rep - optional callback which is called after each iteration | |
| can be NULL | | can be NULL | |
| ptr - optional pointer which is passed to func/grad/hess/jac/rep | | ptr - optional pointer which is passed to func/grad/hess/jac/rep | |
| can be NULL | | can be NULL | |
| | | | |
| NOTES: | | NOTES: | |
| | | | |
| 1. Depending on function used to create state structure, this algorithm | | 1. Depending on function used to create state structure, this algorithm | |
| | | | |
| skipping to change at line 989 | | skipping to change at line 1239 | |
|---|
| will be no callback to call. In this case exception will be thrown. | | will be no callback to call. In this case exception will be thrown. | |
| | | | |
| Be careful to avoid such errors because there is no way to find them at | | Be careful to avoid such errors because there is no way to find them at | |
| compile time - you can see them at runtime only. | | compile time - you can see them at runtime only. | |
| | | | |
| -- ALGLIB -- | | -- ALGLIB -- | |
| Copyright 10.03.2009 by Bochkanov Sergey | | Copyright 10.03.2009 by Bochkanov Sergey | |
| | | | |
| *************************************************************************/ | | *************************************************************************/ | |
| void minlmoptimize(minlmstate &state, | | void minlmoptimize(minlmstate &state, | |
|
| void (*func)(const real_1d_array &x, double &func, void *ptr), | | void (*fvec)(const real_1d_array &x, real_1d_array &fi, void *ptr), | |
| | | void (*rep)(const real_1d_array &x, double func, void *ptr) = NULL, | |
| | | void *ptr = NULL); | |
| | | void minlmoptimize(minlmstate &state, | |
| | | void (*fvec)(const real_1d_array &x, real_1d_array &fi, void *ptr), | |
| void (*jac)(const real_1d_array &x, real_1d_array &fi, real_2d_array &
jac, void *ptr), | | void (*jac)(const real_1d_array &x, real_1d_array &fi, real_2d_array &
jac, void *ptr), | |
| void (*rep)(const real_1d_array &x, double func, void *ptr) = NULL, | | void (*rep)(const real_1d_array &x, double func, void *ptr) = NULL, | |
| void *ptr = NULL); | | void *ptr = NULL); | |
| void minlmoptimize(minlmstate &state, | | void minlmoptimize(minlmstate &state, | |
| void (*func)(const real_1d_array &x, double &func, void *ptr), | | void (*func)(const real_1d_array &x, double &func, void *ptr), | |
| void (*grad)(const real_1d_array &x, double &func, real_1d_array &grad,
void *ptr), | | void (*grad)(const real_1d_array &x, double &func, real_1d_array &grad,
void *ptr), | |
|
| | | void (*hess)(const real_1d_array &x, double &func, real_1d_array &grad, | |
| | | real_2d_array &hess, void *ptr), | |
| | | void (*rep)(const real_1d_array &x, double func, void *ptr) = NULL, | |
| | | void *ptr = NULL); | |
| | | void minlmoptimize(minlmstate &state, | |
| | | void (*func)(const real_1d_array &x, double &func, void *ptr), | |
| void (*jac)(const real_1d_array &x, real_1d_array &fi, real_2d_array &
jac, void *ptr), | | void (*jac)(const real_1d_array &x, real_1d_array &fi, real_2d_array &
jac, void *ptr), | |
| void (*rep)(const real_1d_array &x, double func, void *ptr) = NULL, | | void (*rep)(const real_1d_array &x, double func, void *ptr) = NULL, | |
| void *ptr = NULL); | | void *ptr = NULL); | |
| void minlmoptimize(minlmstate &state, | | void minlmoptimize(minlmstate &state, | |
| void (*func)(const real_1d_array &x, double &func, void *ptr), | | void (*func)(const real_1d_array &x, double &func, void *ptr), | |
| void (*grad)(const real_1d_array &x, double &func, real_1d_array &grad,
void *ptr), | | void (*grad)(const real_1d_array &x, double &func, real_1d_array &grad,
void *ptr), | |
|
| void (*hess)(const real_1d_array &x, double &func, real_1d_array &grad,
real_2d_array &hess, void *ptr), | | void (*jac)(const real_1d_array &x, real_1d_array &fi, real_2d_array &
jac, void *ptr), | |
| void (*rep)(const real_1d_array &x, double func, void *ptr) = NULL, | | void (*rep)(const real_1d_array &x, double func, void *ptr) = NULL, | |
| void *ptr = NULL); | | void *ptr = NULL); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
| Levenberg-Marquardt algorithm results | | Levenberg-Marquardt algorithm results | |
| | | | |
| INPUT PARAMETERS: | | INPUT PARAMETERS: | |
| State - algorithm state | | State - algorithm state | |
| | | | |
| OUTPUT PARAMETERS: | | OUTPUT PARAMETERS: | |
| X - array[0..N-1], solution | | X - array[0..N-1], solution | |
|
| Rep - optimization report: | | Rep - optimization report; | |
| * Rep.TerminationType completetion code: | | see comments for this structure for more info. | |
| * -1 incorrect parameters were specified | | | |
| * 1 relative function improvement is no more than | | | |
| EpsF. | | | |
| * 2 relative step is no more than EpsX. | | | |
| * 4 gradient is no more than EpsG. | | | |
| * 5 MaxIts steps was taken | | | |
| * 7 stopping conditions are too stringent, | | | |
| further improvement is impossible | | | |
| * Rep.IterationsCount contains iterations count | | | |
| * Rep.NFunc - number of function calculations | | | |
| * Rep.NJac - number of Jacobi matrix calculations | | | |
| * Rep.NGrad - number of gradient calculations | | | |
| * Rep.NHess - number of Hessian calculations | | | |
| * Rep.NCholesky - number of Cholesky decomposition calculat | | | |
| ions | | | |
| | | | |
| -- ALGLIB -- | | -- ALGLIB -- | |
| Copyright 10.03.2009 by Bochkanov Sergey | | Copyright 10.03.2009 by Bochkanov Sergey | |
| *************************************************************************/ | | *************************************************************************/ | |
| void minlmresults(const minlmstate &state, real_1d_array &x, minlmreport &r
ep); | | void minlmresults(const minlmstate &state, real_1d_array &x, minlmreport &r
ep); | |
| | | | |
| /************************************************************************* | | /************************************************************************* | |
| Levenberg-Marquardt algorithm results | | Levenberg-Marquardt algorithm results | |
| | | | |
| Buffered implementation of MinLMResults(), which uses pre-allocated buffer | | Buffered implementation of MinLMResults(), which uses pre-allocated buffer | |
| | | | |
| skipping to change at line 1524 | | skipping to change at line 1769 | |
|---|
| ae_state *_state); | | ae_state *_state); | |
| void minlbfgssetstpmax(minlbfgsstate* state, | | void minlbfgssetstpmax(minlbfgsstate* state, | |
| double stpmax, | | double stpmax, | |
| ae_state *_state); | | ae_state *_state); | |
| void minlbfgscreatex(ae_int_t n, | | void minlbfgscreatex(ae_int_t n, | |
| ae_int_t m, | | ae_int_t m, | |
| /* Real */ ae_vector* x, | | /* Real */ ae_vector* x, | |
| ae_int_t flags, | | ae_int_t flags, | |
| minlbfgsstate* state, | | minlbfgsstate* state, | |
| ae_state *_state); | | ae_state *_state); | |
|
| | | void minlbfgssetdefaultpreconditioner(minlbfgsstate* state, | |
| | | ae_state *_state); | |
| | | void minlbfgssetcholeskypreconditioner(minlbfgsstate* state, | |
| | | /* Real */ ae_matrix* p, | |
| | | ae_bool isupper, | |
| | | ae_state *_state); | |
| ae_bool minlbfgsiteration(minlbfgsstate* state, ae_state *_state); | | ae_bool minlbfgsiteration(minlbfgsstate* state, ae_state *_state); | |
| void minlbfgsresults(minlbfgsstate* state, | | void minlbfgsresults(minlbfgsstate* state, | |
| /* Real */ ae_vector* x, | | /* Real */ ae_vector* x, | |
| minlbfgsreport* rep, | | minlbfgsreport* rep, | |
| ae_state *_state); | | ae_state *_state); | |
| void minlbfgsresultsbuf(minlbfgsstate* state, | | void minlbfgsresultsbuf(minlbfgsstate* state, | |
| /* Real */ ae_vector* x, | | /* Real */ ae_vector* x, | |
| minlbfgsreport* rep, | | minlbfgsreport* rep, | |
| ae_state *_state); | | ae_state *_state); | |
| void minlbfgsrestartfrom(minlbfgsstate* state, | | void minlbfgsrestartfrom(minlbfgsstate* state, | |
| /* Real */ ae_vector* x, | | /* Real */ ae_vector* x, | |
| ae_state *_state); | | ae_state *_state); | |
| ae_bool _minlbfgsstate_init(minlbfgsstate* p, ae_state *_state, ae_bool mak
e_automatic); | | ae_bool _minlbfgsstate_init(minlbfgsstate* p, ae_state *_state, ae_bool mak
e_automatic); | |
| ae_bool _minlbfgsstate_init_copy(minlbfgsstate* dst, minlbfgsstate* src, ae
_state *_state, ae_bool make_automatic); | | ae_bool _minlbfgsstate_init_copy(minlbfgsstate* dst, minlbfgsstate* src, ae
_state *_state, ae_bool make_automatic); | |
| void _minlbfgsstate_clear(minlbfgsstate* p); | | void _minlbfgsstate_clear(minlbfgsstate* p); | |
| ae_bool _minlbfgsreport_init(minlbfgsreport* p, ae_state *_state, ae_bool m
ake_automatic); | | ae_bool _minlbfgsreport_init(minlbfgsreport* p, ae_state *_state, ae_bool m
ake_automatic); | |
| ae_bool _minlbfgsreport_init_copy(minlbfgsreport* dst, minlbfgsreport* src,
ae_state *_state, ae_bool make_automatic); | | ae_bool _minlbfgsreport_init_copy(minlbfgsreport* dst, minlbfgsreport* src,
ae_state *_state, ae_bool make_automatic); | |
| void _minlbfgsreport_clear(minlbfgsreport* p); | | void _minlbfgsreport_clear(minlbfgsreport* p); | |
|
| | | void minlmcreatevj(ae_int_t n, | |
| | | ae_int_t m, | |
| | | /* Real */ ae_vector* x, | |
| | | minlmstate* state, | |
| | | ae_state *_state); | |
| | | void minlmcreatev(ae_int_t n, | |
| | | ae_int_t m, | |
| | | /* Real */ ae_vector* x, | |
| | | double diffstep, | |
| | | minlmstate* state, | |
| | | ae_state *_state); | |
| void minlmcreatefgh(ae_int_t n, | | void minlmcreatefgh(ae_int_t n, | |
| /* Real */ ae_vector* x, | | /* Real */ ae_vector* x, | |
| minlmstate* state, | | minlmstate* state, | |
| ae_state *_state); | | ae_state *_state); | |
|
| | | void minlmcreatevgj(ae_int_t n, | |
| | | ae_int_t m, | |
| | | /* Real */ ae_vector* x, | |
| | | minlmstate* state, | |
| | | ae_state *_state); | |
| void minlmcreatefgj(ae_int_t n, | | void minlmcreatefgj(ae_int_t n, | |
| ae_int_t m, | | ae_int_t m, | |
| /* Real */ ae_vector* x, | | /* Real */ ae_vector* x, | |
| minlmstate* state, | | minlmstate* state, | |
| ae_state *_state); | | ae_state *_state); | |
| void minlmcreatefj(ae_int_t n, | | void minlmcreatefj(ae_int_t n, | |
| ae_int_t m, | | ae_int_t m, | |
| /* Real */ ae_vector* x, | | /* Real */ ae_vector* x, | |
| minlmstate* state, | | minlmstate* state, | |
| ae_state *_state); | | ae_state *_state); | |
| void minlmsetcond(minlmstate* state, | | void minlmsetcond(minlmstate* state, | |
| double epsg, | | double epsg, | |
| double epsf, | | double epsf, | |
| double epsx, | | double epsx, | |
| ae_int_t maxits, | | ae_int_t maxits, | |
| ae_state *_state); | | ae_state *_state); | |
| void minlmsetxrep(minlmstate* state, ae_bool needxrep, ae_state *_state); | | void minlmsetxrep(minlmstate* state, ae_bool needxrep, ae_state *_state); | |
| void minlmsetstpmax(minlmstate* state, double stpmax, ae_state *_state); | | void minlmsetstpmax(minlmstate* state, double stpmax, ae_state *_state); | |
|
| | | void minlmsetacctype(minlmstate* state, | |
| | | ae_int_t acctype, | |
| | | ae_state *_state); | |
| ae_bool minlmiteration(minlmstate* state, ae_state *_state); | | ae_bool minlmiteration(minlmstate* state, ae_state *_state); | |
| void minlmresults(minlmstate* state, | | void minlmresults(minlmstate* state, | |
| /* Real */ ae_vector* x, | | /* Real */ ae_vector* x, | |
| minlmreport* rep, | | minlmreport* rep, | |
| ae_state *_state); | | ae_state *_state); | |
| void minlmresultsbuf(minlmstate* state, | | void minlmresultsbuf(minlmstate* state, | |
| /* Real */ ae_vector* x, | | /* Real */ ae_vector* x, | |
| minlmreport* rep, | | minlmreport* rep, | |
| ae_state *_state); | | ae_state *_state); | |
| void minlmrestartfrom(minlmstate* state, | | void minlmrestartfrom(minlmstate* state, | |
| | | | |
| End of changes. 61 change blocks. |
|---|
| 96 lines changed or deleted | | 373 lines changed or added | |
|---|
|