functions.cxx | searchcode

/share/functions.cxx

https://gitlab.com/ineris/amc
C++ | 570 lines | 396 code | 108 blank | 66 comment | 84 complexity | 2011047747ca3434140b72312d66e57b MD5 | raw file
Possible License(s): GPL-3.0

// Copyright (C) 2013-2015 INERIS
// Author(s) : Edouard Debry
// 
// This file is part of AMC.
// 
// AMC is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
// 
// AMC is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
// 
// You should have received a copy of the GNU General Public License
// along with AMC.  If not, see <http://www.gnu.org/licenses/>.

#ifndef AMC_FILE_FUNCTIONS_CXX

#include "functions.hxx"

namespace AMC
{
  // Compute cubic root.
  inline float pow3(float x)
  {
    // Save constant befor modifying it.
    float a = x * 0.3333333333333333333;

    // First approximation.
    const int ebits = 8;
    const int fbits = 23;
    const int n = 3;

    int32_t& i = (int32_t&) x;
    const int32_t bias = (1 << (ebits - 1)) - 1;
    i = (i - (bias << fbits)) / n + (bias << fbits);

    x = 0.6666666666666666666 * x + a / (x * x);
    x = 0.6666666666666666666 * x + a / (x * x);
    x = 0.6666666666666666666 * x + a / (x * x);
    x = 0.6666666666666666666 * x + a / (x * x);

    return x;
  }


  inline double pow3(double x)
  {
    // Save constant befor modifying it.
    double a = x * 0.3333333333333333333;

    // First approximation.
    const int ebits = 11;
    const int fbits = 52;
    const int n = 3;

    int64_t& i = (int64_t&) x;
    const int64_t bias = (1 << (ebits - 1)) - 1;
    i = (i - (bias << fbits)) / n + (bias << fbits);

    x = 0.6666666666666666666 * x + a / (x * x);
    x = 0.6666666666666666666 * x + a / (x * x);
    x = 0.6666666666666666666 * x + a / (x * x);
    x = 0.6666666666666666666 * x + a / (x * x);

    return x;
  }


  // Get the number of power of p out of n.
  inline int get_power_out_of(const int &n, const int p)
  {
    int i(0), j(n);
    while (j / p > 0)
      {
        j /= p;
        i++;
      }

    return i;
  }


  // Check index.
  inline void check_index(const int &index_min, const int &index_max, const int &index)
  {
    if (index < index_min || index >= index_max)
      throw string("Out of range index, should be within ["
                   + to_str(index_min) + ", " + to_str(index_max)
                   + "[, got " + to_str(index) + ".");
  }


  // Get combinations of p values within n.
  bool get_combination(const int &p, const int &n, vector<int> &index_reference, vector<int> &index)
  {
    int i, j, k;

    if (int(index_reference.size()) == 0)
      for (i = 0; i < p; i++)
        index_reference.push_back(i);

    if (index_reference == index)
      return false;

    index = index_reference;

    for (i = 0; i < p; i++)
      {
        j = p - i - 1;
        if (index_reference[j] < n - 1 - i)
          {
            index_reference[j]++;
            for (k = j + 1; k < p; k++)
              index_reference[k] = index_reference[k - 1] + 1;       
            break;
          }
      }

    return true;
  }


  // Find convex hull of a polygon.
  template<class T>
  void find_polygon_convex_hull(const int &p, const int &n, const vector<T> &point, vector<int> &edge)
  {
    // n is the number of points.
    // p is the points dimension.
    // Then, point should be of dimension n * p.
    if (int(point.size()) != (n * p))
      throw string("Size of point vector is " + to_str(point.size()) + ", expected " + to_str(n * p) + ".");

    int q = p + 1;
    vector<int> index_reference, index;

    // When number of points is just below the dimension all points are edges.
    if (n == q)
      while(get_combination(p, n, index_reference, index))
        for (int i = 0; i < p; i++)
          edge.push_back(index[i]);
    else
      {
        int p2 = p * p;
        int q2 = q * q;

        while(get_combination(p, n, index_reference, index))
          {
            // Compute supplementary index.
            vector<int> index_supplementary;
            for (int i = 0; i < n; i++)
              {
                bool is_included(false);
                for (int j = 0; j < p; j++)
                  if (i == index[j])
                    {
                      is_included = true;
                      break;
                    }

                if (! is_included)
                  index_supplementary.push_back(i);
              }

            int np = int(index_supplementary.size());
            if (np != (n - p))
              throw string("Wrong number of supplementary points, got " + to_str(np) + ", expected " + to_str(n - p) + ".");

            // Compute intermediate determinants.
            vector<T> mat(q2, T(1));

            int i(q);
            for (int j = 0; j < p; j++)
              {
                int l = p * index[j];
                for (int k = 0; k < p; k++)
                  mat[i++] = point[l + k];
                i++;
              }

            vector<T> det(np);
            for (i = 0; i < np; i++)
              {
                int l = p * index_supplementary[i];

                for (int j = 0; j < p; j++)
                  mat[j] = point[l + j];

                det[i] = compute_matrix_determinant<T>(mat, q);
              }

            for (i = 0; i < np; i++)
              if (det[i] == T(0))
                throw string("Got coplanar points.");

            bool same_sign(true);
            for (i = 1; i < np; i++)
              if (det[0] * det[i] < T(0))
                {
                  same_sign = false;
                  break;
                }

            // All determinants are of same sign,
            // this means all supplementary points are on
            // the same side of given surface, provided
            // the polygon is convex, this entails the surface
            // to be an edge.
            if (same_sign)
              for (i = 0; i < p; i++)
                edge.push_back(index[i]);
          }
      }
  }


  // Search index in a vector.
  template<typename T>
  inline int search_index(const vector<T> &v, const T &x, int imin, int imax)
  {
    imax = imax < 0 ? int(v.size()) + imax : imax;

    // Search algorithm
    if (x < v[imin])
      return imin;
    else if (x >= v[imax])
      return imax;
    else
      {
        int imid = (imin + imax) / 2;

        while (imax >= imin)
          {
            if (x > v[imid])
              imin = imid + 1;
            else if (x < v[imid])
              imax = imid - 1;
            else
              return imid;

            imid = (imin + imax) / 2;
          }

        return imid;
      }
  }


  template<typename T>
  inline int search_index_ascending(const vector<T> &v, const T x, int imin)
  {
    int imax = int(v.size()) - 1;

    for (int i = ++imin; i < imax; ++i)
      if (x < v[i])
        return --i;

    return --imax;
  }


  template<typename T>
  inline int search_index_descending(const vector<T> &v, const T x, int imax)
  {
    imax = imax < 0 ? int(v.size()) - 2 : imax;

    for (int i = imax; i > 0; --i)
      if (x >= v[i])
        return i;

    return 0;
  }


  // Test is one vector (v) is included into another one (w).
  template<typename T>
  bool is_included(const vector<T> &v, const vector<T> &w)
  {
    if (v.size() > w.size())
      throw string("Vector v has a greater size as vector w, hence cannot be included in this one.");

    size_t count(0);
    for (int i = 0; i < int(v.size()); i++)
      for (int j = 0; j < int(w.size()); j++)
        if (v[i] == w[j])
          {
            count++;
            break;
          }

    return count == v.size();
  }


  // Compute sum of vector or part of vector.
  template<typename T>
  inline T compute_vector_sum(const vector<T> &v, int n)
  {
    if (n <= 0)
      n += v.size();
    
    T sum(T(0));
    for (int i = 0; i < n; i++)
      sum += v[i];
    return sum;
  }


  // Compute the determinant of a matrix, stored in row major order (C-style).
  template<typename T>
  T compute_matrix_determinant(const vector<T> &matrix, int n)
  {
    if (n == 0)
      return 0;
    else if (n == 1)
      return matrix[0];
    else if (n == 2)
      return matrix[0] * matrix[3] - matrix[1] * matrix[2];
    else
      {
        T determinant(T(0));

        for (int i = 0; i < n; i++)
          {
            int n2 = n - 1;
            vector<T> matrix2(n2 * n2);

            int l(0), m(n);
            for (int j = 1; j < n; j++)
              {
                for (int k = 0; k < i; k++)
                  matrix2[l++] = matrix[m++];
                m++;
                for (int k = i + 1; k < n; k++)
                  matrix2[l++] = matrix[m++];
              }

            determinant += matrix[i]
              * ((i % 2 == 0) ? T(1) : T(-1))
              * compute_matrix_determinant(matrix2, n2);
          }

        return determinant;
      }
  }


  // Generate pseudo-random vector with law uniform.
  template<typename T>
  inline void generate_random_vector(const int &n, vector<T> &random, const T sum)
  {
    random.assign(n, sum);
    for (int i = 0; i < n - 1; i++)
      random[i] = sum * T(drand48());
    sort(random.begin(), random.end());

    for (int i = n - 1; i > 0; --i)
      random[i] -= random[i - 1];
  }


  // Set the random generator seed in C++.
  void set_random_generator(const int &seed)
  {
    srand48(seed);
  }


  // Generate one random vector.
  template<typename T>
  vector<T> generate_random_vector(const int &n, const T sum)
  {
    vector<T> random;
    generate_random_vector(n, random, sum);
    return random;
  }


  // Split one 2D domain into most equally possible subdomains.
  void split_domain(const int n, int &p, int &q)
  {
    vector<int> vp, vq, rank;

    const int Np(p), Nq(q);

    p = 0;
    q = 0;
    while (p < n)
          {
            p++;

            if (n % p == 0)
              {
                vp.push_back(p);
                q = n / p;
                vq.push_back(q);
                rank.push_back(abs(Nq * p - Np * q));
              }
          }

    // Which couple is the best ?
    int i(0);
    for (int j = 1; j < int(rank.size()); ++j)
      if (rank[j] < rank[i]) i = j;

    p = vp[i];
    q = vq[i];
  }


  // Compute the Gauss-Legendre quadrature zeros and weights, from Numerical Recipes in C.
  template<typename T>
  void compute_gauss_legendre_zero_weight(const int &n, T *x, T *w,
                                          const T x1, const T x2,
                                          const T eps)
  {
    // High precision is a good idea for this routine.
    int m, j, i;
    double z1, z, xm, xl, pp, p3, p2, p1;

    // The roots are symmetric in the interval, so we only have to find half of them.
    m = (n + 1) / 2;
    xm = 0.5 * (x2 + x1);
    xl = 0.5 * (x2 - x1);

    // Loop over the desired roots.
    for (i = 1; i <= m; i++)
      {
        z = cos(3.14159265358979323846 * (i - 0.25) / (n + 0.5));
        // Starting with the above approximation to the ith root,
        // we enter the main loop of reﬁnement by Newton’s method.
        do
          {
            p1 = 1.0;
            p2 = 0.0;
            // Loop up the recurrence relation to get the Legendre polynomial evaluated at z.
            for (j = 1; j <= n; j++)
              {
                p3 = p2;
                p2 = p1;
                p1 = ((2.0 * j - 1.0) * z * p2 - (j - 1.0) * p3) / j;
              }
            // p1 is now the desired Legendre polynomial. We next compute pp, its derivative,
            // by a standard relation involving also p2, the polynomial of one lower order.
            pp = n * (z * p1 - p2) / (z * z - 1.0);
            z1 = z;
            z = z1 - p1 / pp;
            // Newton’s method.
          }
        while (fabs(z - z1) > eps);

        // Scale the root to the desired interval.
        x[i] = xm - xl * z;
        x[n+1-i] = xm + xl * z;
        w[i] = 2.0 * xl / ((1.0 - z * z) * pp * pp);
        // Compute the weight and its symmetric counterpart.
        w[n + 1 - i] = w[i];
      }
  }


  // Comute cubic splines.
  template<typename T>
  void compute_cubic_spline(const int n, const T *x, const T *y, T *y2, const T *yp1, const T *ypn)
  {
    vector<T> u(n - 1, T(0));

    if (yp1 == NULL)
      {
        y2[0] = T(0);
        u[0] = T(0);
      }
    else
      {
        y2[0] = -T(0.5);
        u[0] = (T(3) / (x[1] - x[0])) * ((y[1] - y[0]) / (x[1] - x[0]) - *yp1);
      }

    for (int i = 1; i < n - 1; ++i)
      {
        const T sig=(x[i] - x[i - 1]) / (x[i + 1] - x[i - 1]);
        const T p = sig * y2[i - 1] + T(2);
        y2[i] = (sig - T(1)) / p;

        u[i] = (y[i + 1] - y[i]) / (x[i + 1] - x[i]) - (y[i] - y[i - 1]) / (x[i] - x[i - 1]);
        u[i]=(T(6) * u[i] / (x[i + 1] - x[i - 1]) - sig * u[i - 1]) / p;
      }

    T qn(T(0)), un(T(0));
    if (ypn != NULL)
      {
        qn = T(0.5);
        un = (T(3) / (x[n - 1] - x[n - 2])) * (*ypn - (y[n - 1] - y[n - 2]) / (x[n - 1] - x[n - 2]));
      }

    y2[n - 1] = (un - qn * u[n - 2]) / (qn * y2[n - 2] + T(1));

    for (int k = n - 2; k >= 0; --k)
      y2[k] = y2[k] * y2[k + 1] + u[k];
  }

  template<typename T>
  T apply_cubic_spline(const int n, const T *xa, const T *ya, const T *y2a, const T &x)
  {
    int klo(0), khi(n - 1);

    while (khi - klo > 1)
      {
        const int k = (khi + klo) >> 1;

        if (xa[k] > x)
          khi = k;
        else
          klo = k;
      }

    const T h = xa[khi] - xa[klo];
    const T a = (xa[khi] - x) / h;
    const T b = (x - xa[klo]) / h;

    return a * ya[klo] + b * ya[khi] + (a * (a * a - T(1)) * y2a[klo] + b * (b * b - T(1)) * y2a[khi]) * (h * h) / T(6);
  }


  // Convert ppm to molecule per cm3.
  template<typename T>
  T convert_ppm_to_molecule_per_cm3(const T &pressure,
                                    const T &temperature,
                                    const T concentration)
  {
    // Na / R / 1.e-6 = 1.e6 / kb.
    return T(7.24297156525251e10) * concentration * pressure / temperature;
  }


  // Convert ppt to molecule per cm3.
  template<typename T>
  T convert_ppt_to_molecule_per_cm3(const T &pressure,
                                    const T &temperature,
                                    const T concentration)
  {
    // Na / R / 1.e-9 = 1.e9 / kb.
    return T(7.24297156525251e7) * concentration * pressure / temperature;
  }


  // Convert molecule per cm3 to ppm.
  template<typename T>
  T convert_molecule_per_cm3_to_ppm(const T &pressure,
                                    const T &temperature,
                                    const T concentration)
  {
    return T(1.3806488e-11) * concentration * temperature / pressure;
  }


  // Convert molecule per cm3 to ppt.
  template<typename T>
  T convert_molecule_per_cm3_to_ppt(const T &pressure,
                                    const T &temperature,
                                    const T concentration)
  {
    return T(1.3806488e-8) * concentration * temperature / pressure;
  }
}

#define AMC_FILE_FUNCTIONS_CXX
#endif