https://github.com/cran/robCompositions
Tip revision: 775c7e1be4c58aaf8adccdd2b92d07aa9cdc265f authored by Matthias Templ on 14 January 2020, 05:10:03 UTC
version 2.2.0
version 2.2.0
Tip revision: 775c7e1
density_estimation.h
#ifndef PROGETTO_DENSITY_ESTIMATION_HPP
#define PROGETTO_DENSITY_ESTIMATION_HPP
#include <Eigen/Dense>
#include <Eigen/Sparse>
#include <iostream>
#include <iterator>
#include <vector>
#include "bspline.h"
#include "sandia_rules.h"
#include "gauss_points_weights.h"
constexpr double tol = 1e-04;
// NOTE: in densityEstimator::solve() there is the following note:
// NOTE: QR should do this automatically without giving any warning (CHECK)
/*! \file
@brief Parmeters container class and methods to compute the solution of the problem
*/
/*!
@brief Parameter class manager
@note These parameters are related to the problem itself, not to each single row.
*/
class parametersManager
{
protected:
/*! Spline degree */
unsigned int k;
/*! order of derivative in penalization term */
unsigned int l;
/*! penalization parameter */
double alpha;
/*! Number of control points */
unsigned int n;
/*! Number of knots - 2 (number of internal ones) */
unsigned int g;
/*! Number of knots including additional ones G = g+k+1 */
unsigned int G;
/*! spline knots */
std::vector<double> knots;
/*! [u,v] support of splines */
double u, v;
/*! control points */
std::vector<double> xcp;
public:
parametersManager
(const unsigned int kk, const unsigned int ll, const double opt_param):
k(kk), l(ll), alpha(opt_param), n(0), g(0), G(0), u(0), v(0) {};
/*!
@brief Read knots from C-array style and store them in the std::vector<double> member
*/
void
readKnots
(const double * inputKnots, const unsigned int & size);
/*!
@brief Read control points from C-array style and store them in the std::vector<double> member
@param cancel If cancel = j, it does not consider j-th column.
Useful for cross-validation. Disabled by default.
*/
void
readXcp
(const double * inputXcp, const unsigned int & size, const int & cancel = -1);
void
set_alpha
(const double & opt_param)
{
alpha = opt_param;
}
unsigned int
get_G() const { return G; };
double
get_u() const { return u; };
double
get_v() const { return v;};
unsigned int
get_k() const{ return k;};
unsigned int
get_n() const{ return n;};
};
/*!
@brief Parameter class manager
@note These are the whole set of parameters
*/
class densityEstimator: public parametersManager
{
private:
/*! Collocation matrix - dimensions nxG */
Eigen::MatrixXd C;
/*! dimensions GxG */
Eigen::MatrixXd M;
Eigen::SparseMatrix<double> S;
/*! dimensions GxG */
Eigen::SparseMatrix<double> DK;
/*! matrix of the problem to solve - dimensions GxG */
Eigen::MatrixXd P;
/*! constant term vector of the problem to solve - dimensions G */
Eigen::VectorXd p;
/*! solution of the problem: c = P^(-)p - dimensions G */
Eigen::VectorXd c;
/*! B-spline coefficients - dimensions G */
Eigen::VectorXd b;
/*! Weights associated to data.
More in general they allow to take in account the various data accuracy */
Eigen::VectorXd weights;
/*! extended vector of knots - with extra ones
dimension: g + 2k + 2 = G + k + 1 */
std::vector<double> lambda;
/*! extended vector of knots of the l-th derivative */
std::vector<double> lambda_der;
void
fill_C
(const std::vector<double>& cp);
void
fill_M
();
void
fill_DK
();
/*! Compute the S_l matrix for the penalization term */
void
fill_S
();
void
set_lambda
(const std::vector<double>& knots);
void
set_lambda_der
(const std::vector<double> & knots);
public:
densityEstimator
(const parametersManager & input): parametersManager(input)
{};
/*!
@brief Compute the matrices of the problem
*/
void
set_matrix
();
void
set_weights
(const Eigen::Block<Eigen::Map<Eigen::Matrix<double, -1, -1>,0, Eigen::Stride<0, 0> >, 1, -1, false> & row);
/*!
@brief Compute the value of the functional to be minimized.
@param ycp Value of data in the control points.
*/
double
eval_J
(const std::vector<double>& ycp);
/*!
@brief Compute the solution of the linear system.
@param bspline Vector of Bspline coefficients computed from the solution of the system.
@param ycp Value of data in the control points.
*/
void
solve
(Eigen::Block<Eigen::Matrix<double, -1, -1>, 1, -1, false> bspline, const std::vector<double>& ycp);
void
print_all
() const;
void
print_sol
() const;
std::vector<double>
get_lambda() const { return lambda; };
};
#endif //PROGETTO_DENSITY_ESTIMATION_HPP