34 #include "gsl/gsl_sf_gamma.h"
42 namespace interpolation {
49 : _pointDim(0), _polyCoeffs() {
53 : _pointDim(pv._pointDim),
54 _polyCoeffs(pv._polyCoeffs.num_row(), pv._polyCoeffs.num_col(), 0.) {
60 : _pointDim(numberOfInputVariables), _polyCoeffs(polynomialCoefficients) {
65 : _pointDim(0), _polyCoeffs() {
81 for(
unsigned int i=0; i<coeff.size(); i++) {
82 int polyOrder = coeff[i].InVariables().size();
83 for(
unsigned int j=0; j<coeff[i].InVariables().size(); j++)
85 if(coeff[i].OutVariable() > valueDim) valueDim = coeff[i].OutVariable();
86 if(polyOrder > maxPolyOrder) maxPolyOrder = polyOrder;
94 for(
unsigned int j=0; j<coeff.size(); j++)
96 int dim = coeff[j].OutVariable();
112 for(
unsigned int i=0; i<
PointDimension(); i++) pointV(i+1) = point[i];
114 for(
unsigned int i=0; i<
ValueDimension(); i++) value[i] = valueV(i+1);
122 for(
unsigned int i=0; i<
ValueDimension(); i++) value(i+1) = answer(i+1);
129 polyVector(i+1) = 1.;
149 check[check_index] = poly_power;
150 nearby_points.push_back(check);
151 if (check_index+1 == check.size())
153 for (
int i = 1; i < int(poly_power); ++i)
155 for (
int i = 0; i <= int(poly_power); ++i) {
156 check[check_index] = i;
164 "PPSolveFactory::GetNearbyPointsSquares",
165 "Point dimension must be > 0"
172 std::vector<std::vector<int> > nearby_points(1, std::vector<int>(point_dim, 0));
174 int this_poly_order = 0;
175 while (
int(nearby_points.size()) < index+1) {
177 this_poly_order += 1;
205 for (
int i = 0; i < pointDimension; ++i)
221 PrintContainer<std::vector<int> >(out, _polyKeyByPower[
_pointDim-1][i], int_separator, str_separator, length, pad_at_start);
225 template <
class Container >
229 std::stringstream strstr1(
"");
230 std::stringstream strstr2(
"");
231 typename Container::const_iterator it1 = container.begin();
232 typename Container::const_iterator it2 = it1;
233 while(it1!=container.end())
236 if(it1 != container.end() && it2 != container.end())
237 strstr1 << (*it1) << T_separator;
238 else strstr1 << (*it1);
242 if(
int(strstr1.str().size()) > length && length > 0) strstr2 << strstr1.str().substr(1, length);
243 else strstr2 << strstr1.str();
245 if(!pad_at_start) out << strstr2.str();
246 for(
int i=0; i<length - int(strstr2.str().size()); i++) out << str_separator;
247 if(pad_at_start) out << strstr2.str();
252 size_t maxPower = *std::max_element(index.begin(), index.end());
static std::vector< int > IndexByPower(int index, int nInputVariables)
void PrintHeader(std::ostream &out, char int_separator='.', char str_separator=' ', int length=14, bool pad_at_start=true) const
unsigned int ValueDimension() const
static std::vector< std::vector< std::vector< int > > > _polyKeyByPower
void SetCoefficients(int pointDim, MMatrix< double > coeff)
MMatrix< double > _polyCoeffs
static void IndexByPowerRecursive(std::vector< int > check, size_t check_index, size_t poly_power, std::vector< std::vector< int > > &nearby_points)
MVector< double > & MakePolyVector(const MVector< double > &point, MVector< double > &polyVector) const
size_t num_col() const
returns number of columns in the matrix
void F(const double *point, double *value) const
static std::vector< std::vector< std::vector< int > > > _polyKeyByVector
static void PrintContainer(std::ostream &out, const Container &container, char T_separator, char str_separator, int length, bool pad_at_start)
SquarePolynomialVector describes a vector of multivariate polynomials.
MMatrix< double > GetCoefficientsAsMatrix() const
static bool _printHeaders
static std::vector< int > IndexByVector(int index, int nInputVariables)
unsigned int PolynomialOrder() const
unsigned int PointDimension() const
static unsigned int NumberOfPolynomialCoefficients(int pointDimension, int order)
std::ostream & operator<<(std::ostream &out, const Mesh::Iterator &it)
size_t num_row() const
returns number of rows in the matrix