math.h

Go to the documentation of this file.

// -----------------------------------------------------------------------------
//
// Copyright (C) 2021 CERN & University of Surrey for the benefit of the
// BioDynaMo collaboration. All Rights Reserved.
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
//
// See the LICENSE file distributed with this work for details.
// See the NOTICE file distributed with this work for additional information
// regarding copyright ownership.
//
// -----------------------------------------------------------------------------
 
#ifndef CORE_UTIL_MATH_H_
#define CORE_UTIL_MATH_H_
 
#include <array>
#include <cmath>
#include <numeric>
#include <vector>
 
#include <TMath.h>
 
#include "core/container/math_array.h"
#include "core/util/log.h"
#include "core/util/random.h"
 
namespace bdm {
 
struct Math {
  static constexpr real_t kPi = static_cast<real_t>(TMath::Pi());
  static constexpr real_t kInfinity = 1e20;
 
  static real_t ToDegree(real_t rad) { return rad * (180 / kPi); }
  static real_t ToRadian(real_t deg) { return deg * (kPi / 180); }
 
  // Helper function that returns distance (L2 norm) between two positions in 3D
  static real_t GetL2Distance(const Real3& pos1, const Real3& pos2) {
    Real3 dist_array;
    dist_array[0] = pos2[0] - pos1[0];
    dist_array[1] = pos2[1] - pos1[1];
    dist_array[2] = pos2[2] - pos1[2];
    return dist_array.Norm();
  }
 
  template <std::size_t N>
  static MathArray<real_t, N> CrossProduct(const MathArray<real_t, N>& a,
                                           const MathArray<real_t, N>& b) {
    MathArray<real_t, N> result;
    result[0] = a[1] * b[2] - a[2] * b[1];
    result[1] = a[2] * b[0] - a[0] * b[2];
    result[2] = a[0] * b[1] - a[1] * b[0];
    return result;
  }
 
  static Real3 Perp3(const Real3& a, real_t random) {
    Real3 vect_perp;
    if (a[0] == 0.0) {
      vect_perp[0] = 1.0;
      vect_perp[1] = 0.0;
      vect_perp[2] = 0.0;
      vect_perp = RotAroundAxis(vect_perp, 6.35 * random, a);
    } else {
      vect_perp[0] = a[1];
      vect_perp[1] = -a[0];
      vect_perp[2] = 0.0;
      vect_perp.Normalize();
      vect_perp = RotAroundAxis(vect_perp, 6.35 * random, a);
    }
    return vect_perp;
  }
 
  static Real3 RotAroundAxis(const Real3& vector, real_t theta,
                             const Real3& axis) {
    auto naxis = axis;
    naxis.Normalize();
 
    auto temp_1 = naxis * (vector * naxis);
    auto temp_2 = (vector - temp_1) * std::cos(-theta);
    auto temp_3 = CrossProduct(vector, naxis) * std::sin(-theta);
 
    return temp_1 + temp_2 + temp_3;
  }
 
  static real_t AngleRadian(const Real3& a, const Real3& b) {
    return std::acos(a * b / (a.Norm() * b.Norm()));
  }
 
  static Real3 ProjectionOnto(const Real3& a, const Real3& b) {
    real_t k = (a * b) / (b * b);
    return b * k;
  }
 
  static real_t MSE(const std::vector<real_t>& v1,
                    const std::vector<real_t>& v2) {
    if (v1.size() != v2.size()) {
      Log::Fatal("Math::MSE", "vectors must have same length");
    }
    real_t error = 0;
    for (size_t i = 0; i < v1.size(); ++i) {
      auto diff = v2[i] - v1[i];
      error += diff * diff;
    }
    return error / v1.size();
  }
};
 
}  // namespace bdm
 
#endif  // CORE_UTIL_MATH_H_