PolynomialFit.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* DO NOT MODIFY THIS HEADER */
3 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
4 /* */
5 /* (c) 2010 Battelle Energy Alliance, LLC */
7 /* */
8 /* Prepared by Battelle Energy Alliance, LLC */
9 /* Under Contract No. DE-AC07-05ID14517 */
10 /* With the U. S. Department of Energy */
11 /* */
12 /* See COPYRIGHT for full restrictions */
13 /****************************************************************/
14
15 #include "PolynomialFit.h"
16
17 // C++ includes
18 #include <fstream>
19
20 extern "C" void FORTRAN_CALL(dgels)(...);
21
23
24 PolynomialFit::PolynomialFit(std::vector<Real> x,
25  std::vector<Real> y,
26  unsigned int order,
27  bool truncate_order)
28  : _x(x), _y(y), _order(order), _truncate_order(truncate_order)
29 {
30  if (_truncate_order) // && (_x.size() / 10) < _order)
31  {
32  if (_x.size() == 1)
33  _order = 0;
34  else
35  {
36  _order = (_x.size() / 10) + 1;
37
38  if (_order > order)
39  _order = order;
40  }
41  }
42  else if (_x.size() < order)
43  throw std::domain_error(
44  "Polynomial Fit requires an order less than the size of the input vector");
45 }
46
47 void
49 {
50  fillMatrix();
52 }
53
54 void
56 {
57  unsigned int num_rows = _x.size();
58  unsigned int num_cols = _order + 1;
59  _matrix.resize(num_rows * num_cols);
60
61  for (unsigned int col = 0; col <= _order; ++col)
62  {
63  for (unsigned int row = 0; row < num_rows; ++row)
64  {
65  Real value = 1;
66  for (unsigned int i = 0; i < col; ++i)
67  value *= _x[row];
68
69  _matrix[(col * num_rows) + row] = value;
70  }
71  }
72 }
73
74 void
76 {
77  char mode = 'N';
78  int num_rows = _x.size();
79  int num_coeff = _order + 1;
80  int num_rhs = 1;
81  int buffer_size = -1;
82  Real opt_buffer_size;
83  Real * buffer;
84  int return_value = 0;
85
86  // Must copy _y because the call to dgels destroys the original values
87  std::vector<Real> rhs = _y;
88
89  FORTRAN_CALL(dgels)
90  (&mode,
91  &num_rows,
92  &num_coeff,
93  &num_rhs,
94  &_matrix[0],
95  &num_rows,
96  &rhs[0],
97  &num_rows,
98  &opt_buffer_size,
99  &buffer_size,
100  &return_value);
101  if (return_value)
102  throw std::runtime_error("Call to Fortran routine 'dgels' returned non-zero exit code");
103
104  buffer_size = (int)opt_buffer_size;
105
106  buffer = new Real[buffer_size];
107  FORTRAN_CALL(dgels)
108  (&mode,
109  &num_rows,
110  &num_coeff,
111  &num_rhs,
112  &_matrix[0],
113  &num_rows,
114  &rhs[0],
115  &num_rows,
116  buffer,
117  &buffer_size,
118  &return_value);
119  delete[] buffer;
120
121  if (return_value)
122  throw std::runtime_error("Call to Fortran routine 'dgels' returned non-zero exit code");
123
124  _coeffs.resize(num_coeff);
125  for (int i = 0; i < num_coeff; ++i)
126  _coeffs[i] = rhs[i];
127 }
128
129 Real
131 {
132  unsigned int size = _coeffs.size();
133  Real value = 0;
134
135  Real curr_x = 1;
136  for (unsigned int i = 0; i < size; ++i)
137  {
138  value += _coeffs[i] * curr_x;
139  curr_x *= x;
140  }
141  return value;
142 }
143
144 void
145 PolynomialFit::dumpSampleFile(std::string base_name,
146  std::string x_label,
147  std::string y_label,
148  float xmin,
149  float xmax,
150  float ymin,
151  float ymax)
152 {
153  std::stringstream filename, filename_pts;
154  const unsigned char fill_character = '0';
155  const unsigned int field_width = 4;
156
157  filename.fill(fill_character);
158  filename << base_name;
159  filename.width(field_width);
160  filename << _file_number << ".plt";
161
162  filename_pts.fill(fill_character);
163  filename_pts << base_name << "_pts";
164  filename_pts.width(field_width);
165  filename_pts << _file_number << ".dat";
166
167  /* First dump the GNUPLOT file with the Least Squares Equations */
168  std::ofstream out(filename.str().c_str());
169  out.precision(15);
170  out.fill(fill_character);
171
172  out << "set terminal postscript color enhanced\n"
173  << "set output \"" << base_name;
174  out.width(field_width);
175  out << _file_number << ".eps\"\n"
176  << "set xlabel \"" << x_label << "\"\n"
177  << "set ylabel \"" << y_label << "\"\n";
178  if (xmin != 0 && xmax != 0)
179  out << "set xrange [" << xmin << ":" << xmax << "]\n";
180  if (ymin != 0 && ymax != 0)
181  out << "set yrange [" << ymin << ":" << ymax << "]\n";
182  out << "set key left top\n"
183  << "f(x)=";
184
185  for (unsigned int i = 0; i < _coeffs.size(); ++i)
186  {
187  if (i)
188  out << "+";
189
190  out << _coeffs[i];
191  for (unsigned int j = 0; j < i; ++j)
192  out << "*x";
193  }
194  out << "\nplot f(x) with lines, '" << filename_pts.str() << "' using 1:2 title \"Points\"\n";
195  out.close();
196
197  libmesh_assert(_x.size() == _y.size());
198
199  out.open(filename_pts.str().c_str());
200  /* Next dump the data points into a seperate file */
201  for (unsigned int i = 0; i < _x.size(); ++i)
202  out << _x[i] << " " << _y[i] << "\n";
203  out << std::endl;
204
205  ++_file_number;
206  out.close();
207 }
208
209 unsigned int
211 {
212  return _x.size();
213 }
214
215 const std::vector<Real> &
217 {
218  return _coeffs;
219 }
std::vector< Real > _matrix
Definition: PolynomialFit.h:96
void FORTRAN_CALL() dgels(...)
const std::vector< Real > & getCoefficients()
Get a const reference to the coefficients of the least squares fit.
static PetscErrorCode Vec x
PolynomialFit(std::vector< Real > X, std::vector< Real > Y, unsigned int order, bool truncate_order=false)
Definition: PolynomialFit.C:24
static int _file_number
void generate()
This function generates the polynomial fit.
Definition: PolynomialFit.C:48
unsigned int getSampleSize()
This function returns the size of the array holding the points, i.e.
std::vector< Real > _x
Definition: PolynomialFit.h:94
bool _truncate_order
Definition: PolynomialFit.h:99
std::vector< Real > _y
Definition: PolynomialFit.h:95
std::vector< Real > _coeffs
Definition: PolynomialFit.h:97
void doLeastSquares()
This function is the wrapper for the LAPACK dgels function and is called by generate.
Definition: PolynomialFit.C:75
void fillMatrix()
This is a helper function that creates the matrix necessary for the Least Squares algorithm...
Definition: PolynomialFit.C:55
Real sample(Real x)
This function will take an independent variable input and will return the dependent variable based on...
void dumpSampleFile(std::string base_name, std::string x_label="X", std::string y_label="Y", float xmin=0, float xmax=0, float ymin=0, float ymax=0)
This function will dump GNUPLOT input files that can be run to show the data points and function fits...
unsigned int _order
Definition: PolynomialFit.h:98