Advertisement
Urbani

Untitled

Jun 28th, 2016
73
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 7.02 KB | None | 0 0
  1. #include <cmath>
  2. #include <cstdio>
  3. #include <vector>
  4. #include <iostream>
  5. #include <algorithm>
  6. #include <limits>
  7. #include <cstdlib>
  8.  
  9. using namespace std;
  10.  
  11.  
  12. class cubic_spline
  13. {
  14. private:
  15.     // Структура, описывающая сплайн на каждом сегменте сетки
  16.     struct spline_tuple
  17.     {
  18.         double a, b, c, d, x;
  19.     };
  20.  
  21.     spline_tuple *splines; // Сплайн
  22.     std::size_t n; // Количество узлов сетки
  23.  
  24.     void free_mem(); // Освобождение памяти
  25.  
  26. public:
  27.     cubic_spline(); //конструктор
  28.     ~cubic_spline(); //деструктор
  29.  
  30.     // Построение сплайна
  31.     // x - узлы сетки, должны быть упорядочены по возрастанию, кратные узлы запрещены
  32.     // y - значения функции в узлах сетки
  33.     // n - количество узлов сетки
  34.     void build_spline(const double *x, const double *y, std::size_t n);
  35.  
  36.     // Вычисление значения интерполированной функции в произвольной точке
  37.     double f(double x) const;
  38. };
  39.  
  40. cubic_spline::cubic_spline() : splines(NULL)
  41. {
  42.  
  43. }
  44.  
  45. cubic_spline::~cubic_spline()
  46. {
  47.     free_mem();
  48. }
  49.  
  50. void cubic_spline::build_spline(const double *x, const double *y, std::size_t n)
  51. {
  52.     free_mem();
  53.  
  54.     this->n = n;
  55.  
  56.     // Инициализация массива сплайнов
  57.     splines = new spline_tuple[n];
  58.     for (std::size_t i = 0; i < n; ++i)
  59.     {
  60.         splines[i].x = x[i];
  61.         splines[i].a = y[i];
  62.     }
  63.     splines[0].c = 0.;
  64.  
  65.     // Решение СЛАУ относительно коэффициентов сплайнов c[i] методом прогонки для трехдиагональных матриц
  66.     // Вычисление прогоночных коэффициентов - прямой ход метода прогонки
  67.     double *alpha = new double[n - 1];
  68.     double *beta = new double[n - 1];
  69.     double A, B, C, F, h_i, h_i1, z;
  70.     alpha[0] = beta[0] = 0.;
  71.     for (std::size_t i = 1; i < n - 1; ++i)
  72.     {
  73.         h_i = x[i] - x[i - 1], h_i1 = x[i + 1] - x[i];
  74.         A = h_i;
  75.         C = 2. * (h_i + h_i1);
  76.         B = h_i1;
  77.         F = 6. * ((y[i + 1] - y[i]) / h_i1 - (y[i] - y[i - 1]) / h_i);
  78.         z = (A * alpha[i - 1] + C);
  79.         alpha[i] = -B / z;
  80.         beta[i] = (F - A * beta[i - 1]) / z;
  81.     }
  82.  
  83.     splines[n - 1].c = (F - A * beta[n - 2]) / (C + A * alpha[n - 2]);
  84.  
  85.     // Нахождение решения - обратный ход метода прогонки
  86.     for (std::size_t i = n - 2; i > 0; --i)
  87.         splines[i].c = alpha[i] * splines[i + 1].c + beta[i];
  88.  
  89.     // Освобождение памяти, занимаемой прогоночными коэффициентами
  90.     delete[] beta;
  91.     delete[] alpha;
  92.  
  93.     // По известным коэффициентам c[i] находим значения b[i] и d[i]
  94.     for (std::size_t i = n - 1; i > 0; --i)
  95.     {
  96.         double h_i = x[i] - x[i - 1];
  97.         splines[i].d = (splines[i].c - splines[i - 1].c) / h_i;
  98.         splines[i].b = h_i * (2. * splines[i].c + splines[i - 1].c) / 6. + (y[i] - y[i - 1]) / h_i;
  99.     }
  100. }
  101.  
  102. double cubic_spline::f(double x) const
  103. {
  104.     if (!splines)
  105.         return std::numeric_limits<double>::quiet_NaN(); // Если сплайны ещё не построены - возвращаем NaN
  106.  
  107.     spline_tuple *s;
  108.     if (x <= splines[0].x) // Если x меньше точки сетки x[0] - пользуемся первым эл-том массива
  109.         s = splines + 1;
  110.     else if (x >= splines[n - 1].x) // Если x больше точки сетки x[n - 1] - пользуемся последним эл-том массива
  111.         s = splines + n - 1;
  112.     else // Иначе x лежит между граничными точками сетки - производим бинарный поиск нужного эл-та массива
  113.     {
  114.         std::size_t i = 0, j = n - 1;
  115.         while (i + 1 < j)
  116.         {
  117.             std::size_t k = i + (j - i) / 2;
  118.             if (x <= splines[k].x)
  119.                 j = k;
  120.             else
  121.                 i = k;
  122.         }
  123.         s = splines + j;
  124.     }
  125.  
  126.     double dx = (x - s->x);
  127.     return s->a + (s->b + (s->c / 2. + s->d * dx / 6.) * dx) * dx; // Вычисляем значение сплайна в заданной точке.
  128. }
  129.  
  130. void cubic_spline::free_mem()
  131. {
  132.     delete[] splines;
  133.     splines = NULL;
  134. }
  135.  
  136. int mounth_to_int(string mounth){
  137.     if(mounth == "January") return 1;
  138.     if(mounth == "February" ) return 2;
  139.     if(mounth == "March" ) return 3;
  140.     if(mounth == "April" ) return 4;
  141.     if(mounth == "May" ) return 5;
  142.     if(mounth == "June" ) return 6;
  143.     if(mounth == "July" ) return 7;
  144.     if(mounth == "August" ) return 8;
  145.     if(mounth == "September" ) return 9;
  146.     if(mounth == "October" ) return 10;
  147.     if(mounth == "November" ) return 11;
  148.     return 12;
  149. }
  150. int main() {
  151.     /* Enter your code here. Read input from STDIN. Print output to STDOUT */  
  152.     int N;
  153.     cin >> N;
  154.     string line;
  155.     vector<vector<string>> splitline(N, vector<string>(4));
  156.     vector<double> mean;
  157.     vector<double> diff;
  158.     vector<double> dates;
  159.     vector<pair<double,int>> miss;
  160.     vector<int> num_miss;
  161.     splitline.resize(N);
  162.     getline(cin, line);
  163.     getline(cin, line);
  164.     for(int i = 0; i < N; i++){
  165.         cin >> splitline[i][0] >> splitline[i][1] >> splitline[i][2] >> splitline[i][3];
  166.        
  167.     }
  168.     for(int i =0; i < N; i++){
  169.         if(splitline[i][2][0] == 'M' || splitline[i][3][0] == 'M'){
  170.             if(splitline[i][2][0] == 'M') miss.push_back(pair<double,int> (12 * atoi(splitline[i][0].c_str()) + mounth_to_int(splitline[i][1]),2));
  171.             else miss.push_back(pair<double,int> (12 * atoi(splitline[i][0].c_str()) + mounth_to_int(splitline[i][1]),3));
  172.             num_miss.push_back(i);
  173.         }
  174.         else {
  175.             mean.push_back((atof(splitline[i][2].c_str()) + atof(splitline[i][3].c_str()))/2);
  176.             diff.push_back(atof(splitline[i][2].c_str()) - atof(splitline[i][3].c_str() ) );
  177.             dates.push_back(12 * atoi(splitline[i][0].c_str()) + mounth_to_int(splitline[i][1]));
  178.         }
  179.     }
  180.     cubic_spline cubic_mean;
  181.     cubic_spline cubic_diff;;
  182.     cubic_mean.build_spline(dates.data(), mean.data(), mean.size());
  183.     cubic_diff.build_spline(dates.data(), diff.data(), diff.size());
  184.     for(int i = 0; i< miss.size(); i++){
  185.         if(miss[i].second == 2) {
  186.             if(splitline[num_miss[i]][3][0] == 'M'){
  187.                 cout << cubic_mean.f(miss[i].first) + cubic_diff.f(miss[i].first)/2 << endl;
  188.             }
  189.             //else cout << cubic_mean.f(miss[i].first) + cubic_diff.f(miss[i].first)/2 << endl;
  190.             else cout << atof(splitline[num_miss[i]][3].c_str()) + cubic_diff.f(miss[i].first) << endl;
  191.             //else cout << 2 * cubic_mean.f(miss[i].first) - atoi(splitline[i][3].c_str()) << endl;
  192.         }
  193.         else  {
  194.             if(splitline[num_miss[i]][2][0] == 'M')
  195.             cout << cubic_mean.f(miss[i].first) - cubic_diff.f(miss[i].first)/2 << endl;
  196.             else cout << atof(splitline[num_miss[i]][2].c_str()) - cubic_diff.f(miss[i].first) << endl;
  197.         }
  198.     }
  199.     return 0;
  200. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement