Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #include <iostream>
- #include <vector>
- // #include <cstdio>
- // #include <stdio.h>
- #include <stdlib.h>
- #include <math.h>
- #include <sstream>
- #include <iomanip>
- using namespace std;
- FILE* pipe;
- class Matrix {
- protected:
- vector<vector<double>> data;
- int n;
- int m;
- public:
- Matrix() {
- n = 0;
- m = 0;
- }
- Matrix (int x, int y) {
- n = x;
- m = y;
- data.assign(n, vector<double>(m));
- }
- Matrix& operator+ (const Matrix& matrix) {
- if (matrix.m != m || matrix.n != n) {
- throw runtime_error("Error: the dimensional problem occurred");
- }
- Matrix *ans = new Matrix(n, m);
- for (int i = 0; i < n; i++) {
- for (int j = 0; j < m; j++) {
- ans->data[i][j] = data[i][j] + matrix.data[i][j];
- }
- }
- return *ans;
- }
- Matrix& operator- (const Matrix& matrix) {
- if (matrix.m != m || matrix.n != n) {
- throw runtime_error("Error: the dimensional problem occurred");
- }
- Matrix *ans = new Matrix(n, m);
- for (int i = 0; i < n; i++) {
- for (int j = 0; j < m; j++) {
- ans->data[i][j] = data[i][j] - matrix.data[i][j];
- }
- }
- return *ans;
- }
- Matrix& operator* (const Matrix& matrix) {
- if (matrix.n != m) {
- throw runtime_error("Error: the dimensional problem occurred");
- }
- Matrix *ans = new Matrix(n, matrix.m);
- for (int i = 0; i < n; i++) {
- for (int j = 0; j < matrix.m; j++) {
- for (int k = 0; k < m; k++) {
- ans->data[i][j] += data[i][k]*matrix.data[k][j];
- }
- }
- }
- return *ans;
- }
- Matrix& transpose() {
- Matrix *ans = new Matrix(m, n);
- for (int i = 0; i < n; i++) {
- for (int j = 0; j < m; j++) {
- ans->data[j][i] = data[i][j];
- }
- }
- return *ans;
- }
- Matrix& operator| (const Matrix& matrix) {
- if (matrix.n != n) {
- throw runtime_error("Error: the dimensional problem occurred");
- }
- Matrix *ans = new Matrix(n, matrix.m + m);
- for (int i = 0; i < n; i++) {
- for (int j = 0; j < m; j++) {
- ans->data[i][j] = data[i][j];
- }
- for (int j = m; j < matrix.m + m; j++) {
- ans->data[i][j] = matrix.data[i][j-m];
- }
- }
- return *ans;
- }
- vector<double>& operator[] (int x) {
- if (x >= n) {
- throw runtime_error("Error: index out of range");
- }
- return data[x];
- }
- friend istream& operator>> (istream &is, Matrix &matrix);
- friend ostream& operator<< (ostream &os, Matrix &matrix);
- double get(int i, int j) {
- if (i >= n || j >= m) {
- throw out_of_range("Out of bound for get");
- }
- return data[i][j];
- }
- bool is_square() {
- return n == m;
- }
- int get_n() {
- return n;
- }
- int get_m() {
- return m;
- }
- void normalize_diagonally() {
- for (int i = 0; i < n; i++) {
- if (data[i][i] != 0) {
- for (int j = i + 1; j < m; j++) {
- data[i][j] /= data[i][i];
- }
- data[i][i] = 1;
- }
- }
- }
- void separate() {
- // for (int i = 0; i < n; i++) {
- // for (int j = 0; j < m - 1; j++) {
- // // if (abs(data[i][j]) < 0.005) {
- // // // // cout << fixed << abs(data[i][j]) << ' ';
- // // } else
- // // // // cout << fixed << data[i][j] << ' ';
- // }
- // // cout << '\n';
- // }
- // for (int j = 0; j < n; j++) {
- // if (abs(data[j][m - 1]) < 0.005) {
- // // cout << fixed << abs(data[j][m - 1]) << ' ';
- // } else
- // // cout << fixed << data[j][m - 1] << ' ';
- // }
- // // cout << '\n';
- }
- void last() {
- // for (int j = 0; j < n; j++) {
- // if (abs(data[j][m - 1]) < 0.005) {
- // // cout << fixed << abs(data[j][m - 1]) << ' ';
- // } else
- // // cout << fixed << data[j][m - 1] << ' ';
- // }
- // // cout << '\n';
- }
- Matrix& inverse();
- };
- class Vector : public Matrix {
- public:
- Vector(vector<double> input) : Matrix(input.size(), 1) {
- for (int i = 0 ; i < input.size(); i++) {
- data[i][0] = input[i];
- }
- }
- };
- class LeastSquareAproxMatrix : public Matrix {
- public:
- LeastSquareAproxMatrix(int m, vector<double> input) : Matrix(input.size(), m + 1) {
- for (int i = 0; i < input.size(); i++) {
- double pow = 1;
- for (int j = 0; j <= m; j++) {
- data[i][j] = pow;
- pow *= input[i];
- }
- }
- }
- };
- class SquareMatrix : public Matrix {
- public:
- SquareMatrix() {
- n = 0;
- m = 0;
- }
- SquareMatrix(int x) : Matrix(x, x) {};
- friend istream& operator>> (istream &is, SquareMatrix &matrix);
- };
- class IdentityMatrix : public SquareMatrix {
- public:
- IdentityMatrix() {
- n = 0;
- m = 0;
- }
- IdentityMatrix(int x) : SquareMatrix(x) {
- for (int i = 0; i < x; i++) {
- data[i][i] = 1;
- }
- }
- };
- class EliminationMatrix : public IdentityMatrix {
- public:
- EliminationMatrix(Matrix& matrix, int j, int i) : IdentityMatrix(matrix.get_n()) {
- data[j - 1][i - 1] = -matrix.get(j - 1, i - 1) / matrix.get(i - 1, i - 1);
- }
- };
- class PermutationMatrix : public IdentityMatrix {
- public:
- PermutationMatrix(int x, int i, int j) : IdentityMatrix(x) {
- data[j - 1][j - 1] = 0;
- data[i - 1][i - 1] = 0;
- data[i - 1][j - 1] = 1;
- data[j - 1][i - 1] = 1;
- }
- };
- istream& operator>> (istream &is, SquareMatrix& matrix) {
- is >> matrix.n;
- matrix.m = matrix.n;
- matrix.data.assign(matrix.n, vector<double>(matrix.m));
- for (int i = 0; i < matrix.n; i++) {
- for (int j = 0; j < matrix.m; j++) {
- is >> matrix.data[i][j];
- }
- }
- return is;
- }
- istream& operator>> (istream &is, Matrix& matrix) {
- is >> matrix.n >> matrix.m;
- matrix.data.assign(matrix.n, vector<double>(matrix.m));
- for (int i = 0; i < matrix.n; i++) {
- for (int j = 0; j < matrix.m; j++) {
- is >> matrix.data[i][j];
- }
- }
- return is;
- }
- ostream& operator<< (ostream &os, Matrix& matrix) {
- for (int i = 0; i < matrix.n; i++) {
- for (int j = 0; j < matrix.m; j++) {
- if (abs(matrix.data[i][j]) < 0.005) {
- os << fixed << abs(matrix.data[i][j]) << ' ';
- } else
- os << fixed << matrix.data[i][j] << ' ';
- }
- os << '\n';
- }
- return os;
- }
- bool gaus_elimination(Matrix& a, int &step) {
- bool sign = 1;
- for (int i = 0; i < a.get_n(); i++) {
- double m = 0;
- int m_i = 0;
- for (int j = i; j < a.get_n(); j++) {
- if (abs(a.get(j, i)) > abs(m)) {
- m = a.get(j, i);
- m_i = j;
- }
- }
- if (m_i != i) {
- // // cout << fixed << "step #" << step << ": permutation" << '\n';
- PermutationMatrix p(a.get_n(), m_i + 1, i + 1);
- a = p * a;
- // a.separate();
- sign = !sign;
- step++;
- }
- for (int j = i + 1; j < a.get_n(); j++) {
- if (a.get(j, i) != 0) {
- // // cout << fixed << "step #" << step << ": elimination" << '\n';
- step++;
- EliminationMatrix e(a, j + 1, i + 1);
- a = e * a;
- // a.separate();
- }
- }
- }
- return sign;
- }
- void gaus_back_elimination(Matrix& a, int &step) {
- for (int i = a.get_n() - 1; i >= 0; i--) {
- if (a.get(i, i) != 0) {
- for (int j = i - 1 ; j >= 0; j--) {
- if (a.get(j, i) != 0) {
- step++;
- EliminationMatrix e(a, j + 1, i + 1);
- a = e * a;
- a.separate();
- }
- }
- }
- }
- }
- Matrix& Matrix::inverse() {
- Matrix* ans = new Matrix(n, m);
- IdentityMatrix I(n);
- Matrix aug = *this | I;
- int step = 1;
- gaus_elimination(aug, step);
- gaus_back_elimination(aug, step);
- aug.normalize_diagonally();
- for (int j = 0; j < aug.get_n(); j++) {
- for (int i = aug.get_n(); i < aug.get_m(); i++) {
- (*ans)[j][i - aug.get_n()] = aug.get(j, i);
- }
- }
- return *ans;
- }
- class ColumnVector : public Matrix {
- public:
- ColumnVector() {
- n = 0;
- m = 0;
- }
- ColumnVector(int x) {
- n = x;
- m = 1;
- data.resize(n, vector<double>(1));
- }
- friend istream& operator>> (istream &is, ColumnVector &matrix);
- friend ostream& operator<< (ostream &os, ColumnVector &matrix);
- };
- istream& operator>> (istream &is, ColumnVector& matrix) {
- is >> matrix.n;
- matrix.m = 1;
- matrix.data.resize(matrix.n,vector<double>(1));
- for (int i = 0; i < matrix.n; i++) {
- is >> matrix.data[i][0];
- }
- return is;
- }
- ostream& operator<< (ostream &os, ColumnVector& matrix) {
- for (int i = 0; i < matrix.n; i++) {
- os << matrix.data[i][0] << ' ';
- }
- os << '\n';
- return os;
- }
- bool is_singular( Matrix a) {
- for (int i = 0; i < a.get_n(); i++) {
- double m = 0;
- int m_i = 0;
- for (int j = i; j < a.get_n(); j++) {
- if (abs(a.get(j, i)) > abs(m)) {
- m = a.get(j, i);
- m_i = j;
- }
- }
- if (m_i != i) {
- PermutationMatrix p(a.get_n(), m_i + 1, i + 1);
- a = p * a;
- }
- for (int j = i + 1; j < a.get_n(); j++) {
- if (a.get(j, i) != 0) {
- EliminationMatrix e(a, j + 1, i + 1);
- a = e * a;
- }
- }
- }
- for (int i = 0; i < a.get_n(); i++) {
- if (a.get(i, i) == 0) {
- return 0;
- }
- }
- return 1;
- }
- inline string to_s(double val) {
- stringstream ss;
- ss << fixed << setprecision(4) << val;
- return ss.str();
- }
- int main() {
- pipe = popen("/usr/bin/gnuplot -persist", "w");
- double min_x = 1000000;
- double max_x = -1000000;
- double min_y = 100000000;
- double max_y = -1000000;
- if (pipe == NULL) {
- return 0;
- }
- int m;
- cin >> m;
- vector<double> data_a, data_b;
- for (int i = 0; i < m; i++) {
- double a, b;
- cin >> a >> b;
- data_a.push_back(a);
- data_b.push_back(b);
- min_x = min(min_x, a);
- min_y = min(min_y, b);
- max_x= max(max_x, a);
- max_y = max(max_y, b);
- }
- int n;
- cin >> n;
- LeastSquareAproxMatrix A(n, data_a);
- cout << A;
- Vector b(data_b);
- cout << b;
- cout.precision(4);
- Matrix c = (A.transpose() * A).inverse() * A.transpose() * b;
- string s = "f(x)=" +to_s(c.get(c.get_n() - 1, 0)) + "*x**" + to_string(c.get_n() - 1);
- for(int i = c.get_n() - 2; i > 0; i--) {
- s += "+(" + to_s(c.get(i, 0)) + "*x**" + to_string(i) + ")";
- }
- s += "+(" + to_s(c.get(0, 0)) + ")";
- fputs(s.c_str(), pipe);
- fprintf(pipe, "\n");
- fputs(("plot ["+ to_string(int(min_x - 1)) + ":" + to_string(int(max_x+ 1)) + "]["+ to_string(int (min_y - 1)) + ":" + to_string(int (max_y + 1)) + "] '-' using 1:2 title 'points' with points pointtype 5, f(x) title 'approximation' with lines").c_str(), pipe);
- fprintf(pipe, "\n");
- for (int i = 0; i < m; i++) {
- fprintf(pipe, "%f\t%f\n", data_a[i], data_b[i]);
- }
- fprintf(pipe, "%s\n", "e");
- cout << s << '\n';
- cout << c;
- pclose(pipe);
- }
Advertisement
Advertisement