auxiliary.cpp
1 #include <iostream>
2 #include <cmath>
3 #include "auxiliary.h"
4 
5 namespace BioLCCC
6 {
7 
8 void fitSpline(const double *x, const double *y, const int n, double * y2)
9 {
10  double a,b,c,d;
11  double * c1 = new double[n-1];
12  double * d1 = new double[n-1];
13  c1[0] = 0.0;
14  d1[0] = 0.0;
15 
16  for (int i = 1; i < n - 1; i++)
17  {
18  a = x[i] - x[i-1];
19  b = 2.0 * (x[i+1] - x[i-1]);
20  c = x[i+1] - x[i];
21  d = 6.0 * ((y[i+1] - y[i]) / (x[i+1] - x[i])
22  - (y[i] - y[i-1]) / (x[i] - x[i-1]));
23  c1[i] = c / (b - c1[i-1] * a);
24  d1[i] = (d - d1[i-1] * a) / (b - c1[i-1] * a);
25  }
26 
27  y2[n-1] = 0.0;
28  for (int i = n - 2; i >= 0; i--)
29  {
30  y2[i] = d1[i] - c1[i] * y2[i+1];
31  }
32 
33  delete[] c1, d1;
34 }
35 
36 double calculateSpline(const double *x, const double *y, const double * y2,
37  const int n, const double x_in)
38 {
39  int j = 0;
40  int j_up = n - 1;
41  while (j_up > j + 1)
42  {
43  if ((x[j] <= x_in) && (x_in <= x[(j + j_up) / 2]))
44  {
45  j_up = (j + j_up) / 2;
46  }
47  else
48  {
49  j = (j + j_up) / 2;
50  }
51  }
52 
53  double dx = x[j+1] - x[j];
54  double a = (x[j+1] - x_in) / dx;
55  double b = (x_in - x[j]) / dx;
56  return a * y[j] + b * y[j+1]
57  + ((a * a * a - a) * y2[j] + (b * b * b - b) * y2[j + 1])
58  * (dx * dx) / 6.0;
59 }
60 
61 double linInterpolate(const double * x, const double * y, const int n,
62  const double x_in)
63 {
64  for (int i=0; i<n-1; ++i)
65  {
66  if ((x[i] <= x_in) && (x_in <= x[i+1]))
67  {
68  return y[i] + (y[i+1] - y[i]) * (x_in - x[i]) / (x[i+1] - x[i]);
69  }
70  }
71  return y[n];
72 }
73 
74 double polInterpolate(const double * x, const double * y, const int n,
75  const double x_in)
76 {
77  double * p = new double[n];
78  for ( int i=0; i<n; i++)
79  {
80  p[i] = y[i];
81  }
82  for (int i=1; i<n; i++)
83  {
84  for (int j=0; j<n-i; j++)
85  {
86  p[j] = (p[j] * (x_in - x[j+i]) + p[j+1] * (x[j] - x_in))
87  / (x[j] - x[j+i]);
88  }
89  }
90  double output = p[0];
91  delete[] p;
92  return output;
93 }
94 
95 double partPolInterpolate(const double * x, const double * y,
96  const int n, const int n_part, const double x_in)
97 {
98  int k = 0;
99  int k_up = n - 1;
100  while (k_up > k + 1)
101  {
102  if ((x[k] <= x_in) && (x_in <= x[(k + k_up) / 2]))
103  {
104  k_up = (k + k_up) / 2;
105  }
106  else
107  {
108  k = (k + k_up) / 2;
109  }
110  }
111 
112  k = ((k - n_part + 1) > 0) ? (k - n_part + 1) : 0;
113  k = (k + 2 * n_part < n ) ? k : n - 2 * n_part;
114 
115  double * p = new double[n_part*2];
116  for (int i=0; i<n_part*2; i++)
117  {
118  p[i] = y[k+i];
119  }
120  for (int i=1; i<n_part*2; i++)
121  {
122  for (int j=0; j<n_part*2-i; j++)
123  {
124  p[j] = (p[j] * (x_in - x[k+j+i]) + p[j+1] * (x[k+j] - x_in))
125  / (x[k+j] - x[k+j+i]);
126  }
127  }
128  double output = p[0];
129  delete[] p;
130  return output;
131 }
132 
133 void solveMatrixEquation(double * m, double * rhs, const int n)
134 {
135  double temp;
136  bool * reduced = new bool[n];
137  for (int i=0; i<n; i++)
138  {
139  reduced[i] = false;
140  }
141 
142  for (int k=0; k<n; k++)
143  {
144  int pivot_i=0;
145  int pivot_j=0;
146  temp = 0.0;
147 
148  for (int i=0; i<n; i++)
149  {
150  if (!reduced[i])
151  {
152  for (int j=0; j<n; j++)
153  {
154  if ((!reduced[j]) && (fabs(m[i * n + j]) >= temp))
155  {
156  pivot_i = i;
157  pivot_j = j;
158  temp = fabs(m[i * n + j]);
159  }
160  }
161  }
162  }
163 
164  reduced[pivot_j] = true;
165 
166  if (pivot_i != pivot_j)
167  {
168  for (int j=0; j<n; j++)
169  {
170  temp = m[pivot_i * n + j];
171  m[pivot_i * n + j] = m[pivot_j * n + j];
172  m[pivot_j * n + j] = temp;
173  }
174 
175  temp = rhs[pivot_i];
176  rhs[pivot_i] = rhs[pivot_j];
177  rhs[pivot_j] = temp;
178  }
179 
180  if (m[pivot_j * n + pivot_j] == 0.0)
181  {
182  throw("The matrix is singular.");
183  }
184 
185  temp = 1.0/ m[pivot_j * n + pivot_j];
186  for (int j=0; j<n; j++)
187  {
188  m[pivot_j * n + j] *= temp;
189  }
190  rhs[pivot_j] *= temp;
191 
192  for (int i=0; i<n; i++)
193  {
194  if (i != pivot_j)
195  {
196  temp = m[i * n + pivot_j];
197  for (int j=0; j<n; j++)
198  {
199  m[i * n + j] -= m[pivot_j * n + j] * temp;
200  }
201  rhs[i] -= rhs[pivot_j] * temp;
202  }
203  }
204  }
205 
206  delete[] reduced;
207 }
208 
209 void fitPolynomial(double * x, double * y, const int n)
210 {
211  double * matrix = new double[n*n];
212  for (int i=0; i<n; i++)
213  {
214  for (int j=0; j<n; j++)
215  {
216  matrix[i*n + j] = pow(x[i], j);
217  }
218  }
219 
220  solveMatrixEquation(matrix, y, n);
221 }
222 
223 double calculatePolynomial(const double * coeffs, const int n, const double x)
224 {
225  double output = 0;
226  for (int i=0; i<n; i++)
227  {
228  output += coeffs[i] * pow(x, i);
229  }
230  return output;
231 }
232 
233 }
void solveMatrixEquation(double *m, double *rhs, const int n)
Solves a linear matrix equation m * x = rhs for square matrix m of size nxn.
Definition: auxiliary.cpp:133
Apart from classes, BioLCCC contains calculation methods and constants.
Definition: auxiliary.h:4
void fitSpline(const double *x, const double *y, const int n, double *y2)
Calculates the second derivatives of a function.
Definition: auxiliary.cpp:8
void fitPolynomial(double *x, double *y, const int n)
Constructs a polynomial whose values at n points x equals y.
Definition: auxiliary.cpp:209
double calculateSpline(const double *x, const double *y, const double *y2, const int n, const double x_in)
Calculates the value of a function using the cubic spline interpolation.
Definition: auxiliary.cpp:36
double linInterpolate(const double *x, const double *y, const int n, const double x_in)
Calculates the value of a function using a piecewise linear interpolation.
Definition: auxiliary.cpp:61
double polInterpolate(const double *x, const double *y, const int n, const double x_in)
Calculates the value of a function using a polynomial interpolation.
Definition: auxiliary.cpp:74
double calculatePolynomial(const double *coeffs, const int n, const double x)
Calculates the value of a polynomial of (n-1)th power at point x.
Definition: auxiliary.cpp:223
double partPolInterpolate(const double *x, const double *y, const int n, const int n_part, const double x_in)
Calculates the value of a function using a partial polynomial interpolation.
Definition: auxiliary.cpp:95