src/core/auxiliary.cpp
00001 #include <iostream>
00002 #include <cmath>
00003 #include "auxiliary.h"
00004 
00005 namespace BioLCCC
00006 {
00007 
00008 void fitSpline(const double *x, const double *y, const int n, double * y2)
00009 {
00010     double a,b,c,d;
00011     double * c1 = new double[n-1];
00012     double * d1 = new double[n-1];
00013     c1[0] = 0.0;
00014     d1[0] = 0.0;
00015     
00016     for (int i = 1; i < n - 1; i++)
00017     {
00018         a = x[i] - x[i-1];
00019         b = 2.0 * (x[i+1] - x[i-1]);
00020         c = x[i+1] - x[i];
00021         d = 6.0 * ((y[i+1] - y[i]) / (x[i+1] - x[i]) 
00022                     - (y[i] - y[i-1]) / (x[i] - x[i-1]));
00023         c1[i] = c / (b - c1[i-1] * a);
00024         d1[i] = (d - d1[i-1] * a) / (b - c1[i-1] * a);
00025     }
00026 
00027     y2[n-1] = 0.0;
00028     for (int i = n - 2; i >= 0; i--)
00029     {
00030         y2[i] = d1[i] - c1[i] * y2[i+1];
00031     }
00032 
00033     delete[] c1, d1;
00034 }
00035 
00036 double calculateSpline(const double *x, const double *y, const double * y2,
00037     const int n, const double x_in)
00038 {
00039     int j = 0;
00040     int j_up = n - 1;
00041     while (j_up > j + 1) 
00042     {
00043         if ((x[j] <= x_in) && (x_in <= x[(j + j_up) / 2]))
00044         {
00045             j_up = (j + j_up) / 2;
00046         }
00047         else
00048         {
00049             j = (j + j_up) / 2;
00050         }
00051     }
00052 
00053     double dx = x[j+1] - x[j];
00054     double a = (x[j+1] - x_in) / dx;
00055     double b = (x_in - x[j]) / dx;
00056     return a * y[j] + b * y[j+1] 
00057         + ((a * a * a - a) * y2[j] + (b * b * b - b) * y2[j + 1])
00058           * (dx * dx) / 6.0;
00059 }
00060 
00061 double linInterpolate(const double * x, const double * y, const int n,
00062                       const double x_in)
00063 {
00064     for (int i=0; i<n-1; ++i)
00065     {
00066         if ((x[i] <= x_in) && (x_in <= x[i+1]))
00067         {
00068             return y[i] + (y[i+1] - y[i]) * (x_in - x[i]) / (x[i+1] - x[i]);
00069         }
00070     }
00071     return y[n];
00072 }
00073 
00074 double polInterpolate(const double * x, const double * y, const int n, 
00075                       const double x_in)
00076 {
00077     double * p = new double[n];
00078     for ( int i=0; i<n; i++)
00079     {
00080         p[i] = y[i];
00081     }
00082     for (int i=1; i<n; i++)
00083     {
00084         for (int j=0; j<n-i; j++)
00085         {
00086             p[j] = (p[j] * (x_in - x[j+i]) + p[j+1] * (x[j] - x_in))
00087                    / (x[j] - x[j+i]);
00088         }
00089     }
00090     double output = p[0];
00091     delete[] p;
00092     return output;
00093 }
00094 
00095 double partPolInterpolate(const double * x, const double * y, 
00096     const int n, const int n_part, const double x_in)
00097 {
00098     int k = 0;
00099     int k_up = n - 1;
00100     while (k_up > k + 1) 
00101     {
00102         if ((x[k] <= x_in) && (x_in <= x[(k + k_up) / 2]))
00103         {
00104             k_up = (k + k_up) / 2;
00105         }
00106         else
00107         {
00108             k = (k + k_up) / 2;
00109         }
00110     }
00111 
00112     k = ((k - n_part + 1) > 0) ? (k - n_part + 1) : 0;
00113     k = (k + 2 * n_part < n ) ? k : n - 2 * n_part;
00114 
00115     double * p = new double[n_part*2];
00116     for (int i=0; i<n_part*2; i++)
00117     {
00118         p[i] = y[k+i];
00119     }
00120     for (int i=1; i<n_part*2; i++)
00121     {
00122         for (int j=0; j<n_part*2-i; j++)
00123         {
00124             p[j] = (p[j] * (x_in - x[k+j+i]) + p[j+1] * (x[k+j] - x_in))
00125                    / (x[k+j] - x[k+j+i]);
00126         }
00127     }
00128     double output = p[0];
00129     delete[] p;
00130     return output;
00131 }
00132 
00133 void solveMatrixEquation(double * m, double * rhs, const int n)
00134 {
00135     double temp;
00136     bool * reduced = new bool[n];
00137     for (int i=0; i<n; i++)
00138     {
00139         reduced[i] = false;
00140     }
00141 
00142     for (int k=0; k<n; k++) 
00143     {
00144         int pivot_i=0;
00145         int pivot_j=0;
00146         temp = 0.0;
00147 
00148         for (int i=0; i<n; i++)
00149         {
00150             if (!reduced[i])
00151             {
00152                 for (int j=0; j<n; j++)
00153                 {
00154                     if ((!reduced[j]) && (fabs(m[i * n + j]) >= temp))
00155                     {
00156                         pivot_i = i;
00157                         pivot_j = j;
00158                         temp = fabs(m[i * n + j]);
00159                     }
00160                 }
00161             }
00162         }
00163 
00164         reduced[pivot_j] = true;
00165 
00166         if (pivot_i != pivot_j)
00167         {
00168             for (int j=0; j<n; j++)
00169             {
00170                 temp = m[pivot_i * n + j];
00171                 m[pivot_i * n + j] = m[pivot_j * n + j];
00172                 m[pivot_j * n + j] = temp;
00173             }
00174 
00175             temp = rhs[pivot_i];
00176             rhs[pivot_i] = rhs[pivot_j];
00177             rhs[pivot_j] = temp;
00178         }
00179 
00180         if (m[pivot_j * n + pivot_j] == 0.0) 
00181         {
00182             throw("The matrix is singular.");
00183         }
00184 
00185         temp = 1.0/ m[pivot_j * n + pivot_j];
00186         for (int j=0; j<n; j++)
00187         {
00188             m[pivot_j * n + j] *= temp;
00189         }
00190         rhs[pivot_j] *= temp;
00191 
00192         for (int i=0; i<n; i++)
00193         {
00194             if (i != pivot_j)
00195             {
00196                 temp = m[i * n + pivot_j];
00197                 for (int j=0; j<n; j++) 
00198                 {
00199                     m[i * n + j] -= m[pivot_j * n + j] * temp;
00200                 }
00201                 rhs[i] -= rhs[pivot_j] * temp;
00202             }
00203         }
00204     }
00205 
00206     delete[] reduced;
00207 }
00208 
00209 void fitPolynomial(double * x, double * y, const int n) 
00210 {
00211     double * matrix = new double[n*n];
00212     for (int i=0; i<n; i++)
00213     {
00214         for (int j=0; j<n; j++)
00215         {
00216             matrix[i*n + j] = pow(x[i], j);
00217         }
00218     }
00219 
00220     solveMatrixEquation(matrix, y, n);
00221 }    
00222 
00223 double calculatePolynomial(const double * coeffs, const int n, const double x)
00224 {
00225     double output = 0;
00226     for (int i=0; i<n; i++)
00227     {
00228         output += coeffs[i] * pow(x, i);
00229     }
00230     return output;
00231 }
00232 
00233 }
 All Classes Namespaces Functions Variables Enumerations Enumerator