From f3ada0ee814957115e26b23ef656f79b41928c12 Mon Sep 17 00:00:00 2001
From: "Joachim Wuttke (o)" <j.wuttke@fz-juelich.de>
Date: Thu, 28 Jan 2016 13:56:42 +0100
Subject: [PATCH] Remove dead code (eigen2d matrix functions); sort functions
in MathFunctions; rename 'value' -> 'x' or 'z'.
---
Core/Tools/inc/MathFunctions.h | 105 ++---
Core/Tools/src/MathFunctions.cpp | 383 ++++++++----------
.../UnitTests/TestCore/MatrixFunctionsTest.h | 77 ----
.../UnitTests/TestCore/SpecialFunctionsTest.h | 42 +-
Tests/UnitTests/TestCore/main.cpp | 2 -
5 files changed, 253 insertions(+), 356 deletions(-)
delete mode 100644 Tests/UnitTests/TestCore/MatrixFunctionsTest.h
diff --git a/Core/Tools/inc/MathFunctions.h b/Core/Tools/inc/MathFunctions.h
index e60aeaea6c0..8048c90f721 100644
--- a/Core/Tools/inc/MathFunctions.h
+++ b/Core/Tools/inc/MathFunctions.h
@@ -3,8 +3,7 @@
// BornAgain: simulate and fit scattering at grazing incidence
//
//! @file Tools/inc/MathFunctions.h
-//! @brief Define many functions in namespace MathFunctions,
-//! and provide inline implementation for most of them
+//! @brief Defines functions in namespace MathFunctions.
//!
//! @homepage http://www.bornagainproject.org
//! @license GNU General Public License v3 or higher (see COPYING)
@@ -25,81 +24,94 @@
#include <vector>
#include <cmath>
-#include "EigenCore.h"
-
//! Various mathematical functions.
namespace MathFunctions
{
- BA_CORE_API_ double Gaussian(double value, double average, double std_dev);
+// ************************************************************************** //
+// Various functions
+// ************************************************************************** //
- BA_CORE_API_ double IntegratedGaussian(double value, double average, double std_dev);
+ BA_CORE_API_ double StandardNormal(double x);
+ BA_CORE_API_ double Gaussian(double x, double average, double std_dev);
+ BA_CORE_API_ double IntegratedGaussian(double x, double average, double std_dev);
- BA_CORE_API_ double GenerateNormalRandom(double average, double std_dev);
+//! Sine integral function: \f$Si(x)\equiv\int_0^x du \sin(u)/u\f$
+ BA_CORE_API_ double Si(double x);
- BA_CORE_API_ double StandardNormal(double value);
+//! sinc function: \f$sinc(x)\equiv\sin(x)/x\f$
+ BA_CORE_API_ double sinc(double x);
- BA_CORE_API_ double GenerateStandardNormalRandom();
+//! Complex sinc function: \f$sinc(x)\equiv\sin(x)/x\f$
+ BA_CORE_API_ complex_t sinc(const complex_t &z);
- BA_CORE_API_ double GenerateUniformRandom();
+//! Complex tanhc function: \f$tanhc(x)\equiv\tanh(x)/x\f$
+ BA_CORE_API_ complex_t tanhc(const complex_t &z);
+
+ BA_CORE_API_ complex_t Laue(const complex_t &z, size_t N);
+
+//! Calculates the geometric series of z to order N
+ complex_t geometricSum(complex_t z, int exponent);
+
+
+// ************************************************************************** //
+// Bessel functions
+// ************************************************************************** //
//! Bessel function of the first kind and order 0
- BA_CORE_API_ double Bessel_J0(double value);
+ BA_CORE_API_ double Bessel_J0(double x);
//! Bessel function of the first kind and order 1
- BA_CORE_API_ double Bessel_J1(double value);
+ BA_CORE_API_ double Bessel_J1(double x);
//! Bessel function Bessel_J1(x)/x
- BA_CORE_API_ double Bessel_C1(double value);
+ BA_CORE_API_ double Bessel_C1(double x);
//! Complex Bessel function of the first kind and order 0
- BA_CORE_API_ complex_t Bessel_J0(const complex_t &value);
+ BA_CORE_API_ complex_t Bessel_J0(const complex_t &z);
//! Complex Bessel function of the first kind and order 1
- BA_CORE_API_ complex_t Bessel_J1(const complex_t &value);
+ BA_CORE_API_ complex_t Bessel_J1(const complex_t &z);
//! Complex Bessel function Bessel_J1(x)/x
- BA_CORE_API_ complex_t Bessel_C1(const complex_t &value);
-
-//! Sine integral function: \f$Si(x)\equiv\int_0^x du \sin(u)/u\f$
- BA_CORE_API_ double Si(double value);
+ BA_CORE_API_ complex_t Bessel_C1(const complex_t &z);
-//! sinc function: \f$sinc(x)\equiv\sin(x)/x\f$
- BA_CORE_API_ double sinc(double value);
+//! Computes complex Bessel function J0(z), using power series and asymptotic expansion
+ BA_CORE_API_ complex_t Bessel_J0_PowSer(const complex_t &z);
-//! Complex sinc function: \f$sinc(x)\equiv\sin(x)/x\f$
- BA_CORE_API_ complex_t sinc(const complex_t &value);
+//! Computes complex Bessel function J0(z), using power series and asymptotic expansion
+ BA_CORE_API_ complex_t Bessel_J1_PowSer(const complex_t &z);
-//! Complex tanhc function: \f$tanhc(x)\equiv\tanh(x)/x\f$
- BA_CORE_API_ complex_t tanhc(const complex_t &value);
- BA_CORE_API_ complex_t Laue(const complex_t &value, size_t N);
+// ************************************************************************** //
+// Fourier transform and convolution
+// ************************************************************************** //
enum EFFTDirection { FORWARD_FFT, BACKWARD_FFT };
- BA_CORE_API_ std::vector<complex_t > FastFourierTransform(const std::vector<complex_t > &data, EFFTDirection tcase);
- BA_CORE_API_ std::vector<complex_t > FastFourierTransform(const std::vector<double > &data, EFFTDirection tcase);
+ BA_CORE_API_ std::vector<complex_t >
+ FastFourierTransform(const std::vector<complex_t > &data, EFFTDirection tcase);
- BA_CORE_API_ std::vector<complex_t> ConvolveFFT(const std::vector<double> &signal, const std::vector<double> &resfunc);
+ BA_CORE_API_ std::vector<complex_t >
+ FastFourierTransform(const std::vector<double > &data, EFFTDirection tcase);
-#ifndef GCCXML_SKIP_THIS
-//! computes the norm element-wise
- BA_CORE_API_ Eigen::Matrix2d Norm(const Eigen::Matrix2cd &M);
+ BA_CORE_API_ std::vector<complex_t>
+ ConvolveFFT(const std::vector<double> &signal, const std::vector<double> &resfunc);
-//! computes the absolute value element-wise
- BA_CORE_API_ Eigen::Matrix2d Abs(const Eigen::Matrix2cd &M);
-//! computes the complex conjugate element-wise
- BA_CORE_API_ Eigen::Matrix2cd Conj(const Eigen::Matrix2cd &M);
+// ************************************************************************** //
+// Random number generators
+// ************************************************************************** //
-//! element-wise product
- BA_CORE_API_ Eigen::Matrix2cd ProductByElement(const Eigen::Matrix2cd &left,
- const Eigen::Matrix2cd &right);
+ BA_CORE_API_ double GenerateUniformRandom();
+ BA_CORE_API_ double GenerateStandardNormalRandom();
+ BA_CORE_API_ double GenerateNormalRandom(double average, double std_dev);
-//! take element-wise real value
- BA_CORE_API_ Eigen::Matrix2d Real(const Eigen::Matrix2cd &M);
-#endif
+
+// ************************************************************************** //
+// Workarounds for std functions
+// ************************************************************************** //
BA_CORE_API_ inline bool isnan(double x)
{
@@ -119,15 +131,6 @@ namespace MathFunctions
#endif
}
-//! Computes complex Bessel function J0(z), using power series and asymptotic expansion
- BA_CORE_API_ complex_t Bessel_J0_PowSer(const complex_t &value);
-
-//! Computes complex Bessel function J0(z), using power series and asymptotic expansion
- BA_CORE_API_ complex_t Bessel_J1_PowSer(const complex_t &value);
-
-//! Calculates the geometric series of z to order N
- complex_t geometricSum(complex_t z, int exponent);
-
} // Namespace MathFunctions
diff --git a/Core/Tools/src/MathFunctions.cpp b/Core/Tools/src/MathFunctions.cpp
index ee6049a4a9b..01ab34534d3 100644
--- a/Core/Tools/src/MathFunctions.cpp
+++ b/Core/Tools/src/MathFunctions.cpp
@@ -27,258 +27,126 @@
#include "gsl/gsl_integration.h"
-double MathFunctions::GenerateNormalRandom(double average, double std_dev)
-{
- return GenerateStandardNormalRandom()*std_dev + average;
-}
-
-double MathFunctions::GenerateUniformRandom()
-{
- int random_int = std::rand();
- return (double)random_int / RAND_MAX;
-}
-
-double MathFunctions::Bessel_J0(double value)
-{
- return gsl_sf_bessel_J0(value);
-}
-
-double MathFunctions::Bessel_J1(double value)
-{
- return gsl_sf_bessel_J1(value);
-}
-
-double MathFunctions::Bessel_C1(double value)
-{
- return value > Numeric::double_epsilon ? gsl_sf_bessel_J1(value)/value : 0.5;
-}
+// ************************************************************************** //
+// Various functions
+// ************************************************************************** //
-complex_t MathFunctions::Bessel_J0(const complex_t &value)
+double MathFunctions::StandardNormal(double x)
{
- if(std::imag(value) < Numeric::double_epsilon) {
- return complex_t(Bessel_J0(std::real(value)), 0.0);
- } else {
- return Bessel_J0_PowSer(value);
- }
+ return std::exp(-x * x / 2.0) / std::sqrt(Units::PI2);
}
-complex_t MathFunctions::Bessel_J1(const complex_t &value)
+double MathFunctions::Gaussian(double x, double average, double std_dev)
{
- if(std::imag(value) < Numeric::double_epsilon) {
- return complex_t(Bessel_J1(std::real(value)), 0.0);
- } else {
- return Bessel_J1_PowSer(value);
- }
+ return StandardNormal((x - average) / std_dev) / std_dev;
}
-complex_t MathFunctions::Bessel_C1(const complex_t &value)
+double MathFunctions::IntegratedGaussian(double x, double average, double std_dev)
{
- if(std::imag(value) < Numeric::double_epsilon) {
- double xv = std::real(value);
- return std::abs(xv) > Numeric::double_epsilon ? MathFunctions::Bessel_J1(xv)/xv : 0.5;
- } else {
- return std::abs(value) > Numeric::double_epsilon ? MathFunctions::Bessel_J1(value)/value : 0.5;
- }
+ double normalized_x = (x - average) / std_dev;
+ double root2 = std::sqrt(2.0);
+ return (gsl_sf_erf(normalized_x / root2) + 1.0) / 2.0;
}
-double MathFunctions::Si(double value) // int_0^x du Sin(u)/u
+double MathFunctions::Si(double x) // int_0^x du Sin(u)/u
{
- return gsl_sf_Si(value);
+ return gsl_sf_Si(x);
}
-double MathFunctions::sinc(double value) // Sin(x)/x
+double MathFunctions::sinc(double x) // Sin(x)/x
{
- return gsl_sf_sinc(value/Units::PI);
+ return gsl_sf_sinc(x/Units::PI);
}
-complex_t MathFunctions::sinc(const complex_t &value) // Sin(x)/x
+complex_t MathFunctions::sinc(const complex_t &z) // Sin(x)/x
{
// This is an exception from the rule that we must not test floating-point numbers for equality.
- // For small non-zero values, sin(z) returns quite accurately z or z-z^3/6.
+ // For small non-zero arguments, sin(z) returns quite accurately z or z-z^3/6.
// There is no loss of precision in computing sin(z)/z.
// Therefore there is no need for an expensive test like abs(z)<eps.
- if( value==complex_t(0.,0.) )
+ if( z==complex_t(0.,0.) )
return complex_t(1.0, 0.0);
- return std::sin(value)/value;
+ return std::sin(z)/z;
}
-complex_t MathFunctions::tanhc(const complex_t &value) // tanh(x)/x
+complex_t MathFunctions::tanhc(const complex_t &z) // tanh(x)/x
{
- if(std::abs(value)<Numeric::double_epsilon)
+ if(std::abs(z)<Numeric::double_epsilon)
return complex_t(1.0, 0.0);
- return std::tanh(value)/value;
+ return std::tanh(z)/z;
}
-complex_t MathFunctions::Laue(const complex_t &value, size_t N) // Exp(iNx/2)*Sin((N+1)x)/Sin(x)
+complex_t MathFunctions::Laue(const complex_t &z, size_t N) // Exp(iNx/2)*Sin((N+1)x)/Sin(x)
{
if (N==0)
return complex_t(1.0, 0.0);
- if(std::abs(value)<Numeric::double_epsilon)
+ if(std::abs(z)<Numeric::double_epsilon)
return complex_t(N+1.0, 0.0);
- return std::exp(complex_t(0.0, 1.0)*value*(double)N/2.0)*std::sin(value*(N+1.0)/2.0)/std::sin(value/2.0);
-}
-
-#ifndef GCCXML_SKIP_THIS
-Eigen::Matrix2d MathFunctions::Norm(const Eigen::Matrix2cd &M) {
- Eigen::Matrix2d result;
- result(0,0) = std::norm((complex_t)M(0,0));
- result(0,1) = std::norm((complex_t)M(0,1));
- result(1,0) = std::norm((complex_t)M(1,0));
- result(1,1) = std::norm((complex_t)M(1,1));
- return result;
+ return std::exp(complex_t(0.0, 1.0)*z*(double)N/2.0)*std::sin(z*(N+1.0)/2.0)/std::sin(z/2.0);
}
-Eigen::Matrix2d MathFunctions::Abs(const Eigen::Matrix2cd &M) {
- Eigen::Matrix2d result;
- result(0,0) = std::abs((complex_t)M(0,0));
- result(0,1) = std::abs((complex_t)M(0,1));
- result(1,0) = std::abs((complex_t)M(1,0));
- result(1,1) = std::abs((complex_t)M(1,1));
- return result;
-}
-
-Eigen::Matrix2cd MathFunctions::Conj(const Eigen::Matrix2cd &M) {
- Eigen::Matrix2cd result;
- result(0,0) = std::conj((complex_t)M(0,0));
- result(0,1) = std::conj((complex_t)M(0,1));
- result(1,0) = std::conj((complex_t)M(1,0));
- result(1,1) = std::conj((complex_t)M(1,1));
+complex_t MathFunctions::geometricSum(complex_t z, int exponent)
+{
+ if (exponent < 1) {
+ throw LogicErrorException("MathFunctions::geometricSeries:"
+ " exponent should be > 0");
+ }
+ complex_t result(0.0, 0.0);
+ double nd = (double)exponent;
+ --exponent;
+ while (exponent > 0) {
+ result += std::pow(z, exponent) * (nd - exponent);
+ --exponent;
+ }
return result;
}
-Eigen::Matrix2cd MathFunctions::ProductByElement(
- const Eigen::Matrix2cd &left, const Eigen::Matrix2cd &right) {
- Eigen::Matrix2cd result;
- result(0,0) = left(0,0) * right(0,0);
- result(0,1) = left(0,1) * right(0,1);
- result(1,0) = left(1,0) * right(1,0);
- result(1,1) = left(1,1) * right(1,1);
- return result;
-}
-Eigen::Matrix2d MathFunctions::Real(const Eigen::Matrix2cd &M) {
- Eigen::Matrix2d result;
- result(0,0) = ((complex_t)M(0,0)).real();
- result(0,1) = ((complex_t)M(0,1)).real();
- result(1,0) = ((complex_t)M(1,0)).real();
- result(1,1) = ((complex_t)M(1,1)).real();
- return result;
-}
-#endif
+// ************************************************************************** //
+// Bessel functions
+// ************************************************************************** //
-double MathFunctions::Gaussian(double value, double average, double std_dev)
+double MathFunctions::Bessel_J0(double x)
{
- return StandardNormal((value - average) / std_dev) / std_dev;
+ return gsl_sf_bessel_J0(x);
}
-double MathFunctions::IntegratedGaussian(double value, double average, double std_dev)
+double MathFunctions::Bessel_J1(double x)
{
- double normalized_value = (value - average) / std_dev;
- double root2 = std::sqrt(2.0);
- return (gsl_sf_erf(normalized_value / root2) + 1.0) / 2.0;
+ return gsl_sf_bessel_J1(x);
}
-double MathFunctions::StandardNormal(double value)
+double MathFunctions::Bessel_C1(double x)
{
- return std::exp(-value * value / 2.0) / std::sqrt(Units::PI2);
+ return x > Numeric::double_epsilon ? gsl_sf_bessel_J1(x)/x : 0.5;
}
-double MathFunctions::GenerateStandardNormalRandom() // using GSL
+complex_t MathFunctions::Bessel_J0(const complex_t &z)
{
- gsl_rng *r;
- r = gsl_rng_alloc(gsl_rng_ranlxs2);
- double result = gsl_ran_ugaussian(r);
- gsl_rng_free(r);
- return result;
-}
-
-//! @brief simple (and unoptimized) wrapper function
-//! for the discrete fast Fourier transformation library (fftw3)
-
-std::vector<complex_t> MathFunctions::FastFourierTransform(const std::vector<complex_t> &data,
- MathFunctions::EFFTDirection ftCase)
-{
- double scale(1.);
- size_t npx = data.size();
-
- fftw_complex *ftData = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * npx);
- fftw_complex *ftResult = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * npx);
- memset(ftData, 0, sizeof(fftw_complex) * npx);
- memset(ftResult, 0, sizeof(fftw_complex) * npx);
-
- for (size_t i = 0; i < npx; i++) {
- ftData[i][0] = data[i].real();
- ftData[i][1] = data[i].imag();
- }
-
- fftw_plan plan;
- switch (ftCase) {
- case MathFunctions::FORWARD_FFT:
- plan = fftw_plan_dft_1d((int)npx, ftData, ftResult, FFTW_FORWARD, FFTW_ESTIMATE);
- break;
- case MathFunctions::BACKWARD_FFT:
- plan = fftw_plan_dft_1d((int)npx, ftData, ftResult, FFTW_BACKWARD, FFTW_ESTIMATE);
- scale = 1. / double(npx);
- break;
- default:
- throw std::runtime_error(
- "MathFunctions::FastFourierTransform -> Panic! Unknown transform case.");
- }
-
- fftw_execute(plan);
-
- // saving data for user
- std::vector<complex_t> outData;
- outData.resize(npx);
- for (size_t i = 0; i < npx; i++) {
- outData[i] = scale * complex_t(ftResult[i][0], ftResult[i][1]);
+ if(std::imag(z) < Numeric::double_epsilon) {
+ return complex_t(Bessel_J0(std::real(z)), 0.0);
+ } else {
+ return Bessel_J0_PowSer(z);
}
-
- fftw_destroy_plan(plan);
- fftw_free(ftData);
- fftw_free(ftResult);
-
- return outData;
}
-//! @brief simple (and unoptimized) wrapper function
-//! for the discrete fast Fourier transformation library (fftw3);
-//! transforms real to complex
-
-std::vector<complex_t> MathFunctions::FastFourierTransform(const std::vector<double> &data,
- MathFunctions::EFFTDirection ftCase)
+complex_t MathFunctions::Bessel_J1(const complex_t &z)
{
- std::vector<complex_t> cdata;
- cdata.resize(data.size());
- for (size_t i = 0; i < data.size(); i++) {
- cdata[i] = complex_t(data[i], 0);
+ if(std::imag(z) < Numeric::double_epsilon) {
+ return complex_t(Bessel_J1(std::real(z)), 0.0);
+ } else {
+ return Bessel_J1_PowSer(z);
}
- return MathFunctions::FastFourierTransform(cdata, ftCase);
}
-//! convolution of two real vectors of equal size
-
-std::vector<complex_t> MathFunctions::ConvolveFFT(const std::vector<double> &signal,
- const std::vector<double> &resfunc)
+complex_t MathFunctions::Bessel_C1(const complex_t &z)
{
- if (signal.size() != resfunc.size()) {
- throw std::runtime_error("MathFunctions::ConvolveFFT() -> This convolution works only for "
- "two vectors of equal size. Use Convolve class instead.");
- }
- std::vector<complex_t> fft_signal
- = MathFunctions::FastFourierTransform(signal, MathFunctions::FORWARD_FFT);
- std::vector<complex_t> fft_resfunc
- = MathFunctions::FastFourierTransform(resfunc, MathFunctions::FORWARD_FFT);
-
- std::vector<complex_t> fft_prod;
- fft_prod.resize(fft_signal.size());
- for (size_t i = 0; i < fft_signal.size(); i++) {
- fft_prod[i] = fft_signal[i] * fft_resfunc[i];
+ if(std::imag(z) < Numeric::double_epsilon) {
+ double xv = std::real(z);
+ return std::abs(xv) > Numeric::double_epsilon ? MathFunctions::Bessel_J1(xv)/xv : 0.5;
+ } else {
+ return std::abs(z) > Numeric::double_epsilon ? MathFunctions::Bessel_J1(z)/z : 0.5;
}
-
- std::vector<complex_t> result
- = MathFunctions::FastFourierTransform(fft_prod, MathFunctions::BACKWARD_FFT);
- return result;
}
//! Computes the complex Bessel function J0(z), using standard power series and asymptotic expansion.
@@ -412,18 +280,123 @@ complex_t MathFunctions::Bessel_J1_PowSer(const complex_t &z)
return cj1;
}
-complex_t MathFunctions::geometricSum(complex_t z, int exponent)
+// ************************************************************************** //
+// Fourier transform and convolution
+// ************************************************************************** //
+
+//! @brief simple (and unoptimized) wrapper function
+//! for the discrete fast Fourier transformation library (fftw3)
+
+std::vector<complex_t>
+MathFunctions::FastFourierTransform(const std::vector<complex_t> &data,
+ MathFunctions::EFFTDirection ftCase)
{
- if (exponent < 1) {
- throw LogicErrorException("MathFunctions::geometricSeries:"
- " exponent should be > 0");
+ double scale(1.);
+ size_t npx = data.size();
+
+ fftw_complex *ftData = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * npx);
+ fftw_complex *ftResult = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * npx);
+ memset(ftData, 0, sizeof(fftw_complex) * npx);
+ memset(ftResult, 0, sizeof(fftw_complex) * npx);
+
+ for (size_t i = 0; i < npx; i++) {
+ ftData[i][0] = data[i].real();
+ ftData[i][1] = data[i].imag();
}
- complex_t result(0.0, 0.0);
- double nd = (double)exponent;
- --exponent;
- while (exponent > 0) {
- result += std::pow(z, exponent) * (nd - exponent);
- --exponent;
+
+ fftw_plan plan;
+ switch (ftCase) {
+ case MathFunctions::FORWARD_FFT:
+ plan = fftw_plan_dft_1d((int)npx, ftData, ftResult, FFTW_FORWARD, FFTW_ESTIMATE);
+ break;
+ case MathFunctions::BACKWARD_FFT:
+ plan = fftw_plan_dft_1d((int)npx, ftData, ftResult, FFTW_BACKWARD, FFTW_ESTIMATE);
+ scale = 1. / double(npx);
+ break;
+ default:
+ throw std::runtime_error(
+ "MathFunctions::FastFourierTransform -> Panic! Unknown transform case.");
+ }
+
+ fftw_execute(plan);
+
+ // saving data for user
+ std::vector<complex_t> outData;
+ outData.resize(npx);
+ for (size_t i = 0; i < npx; i++) {
+ outData[i] = scale * complex_t(ftResult[i][0], ftResult[i][1]);
+ }
+
+ fftw_destroy_plan(plan);
+ fftw_free(ftData);
+ fftw_free(ftResult);
+
+ return outData;
+}
+
+//! @brief simple (and unoptimized) wrapper function
+//! for the discrete fast Fourier transformation library (fftw3);
+//! transforms real to complex
+
+std::vector<complex_t>
+MathFunctions::FastFourierTransform(const std::vector<double> &data,
+ MathFunctions::EFFTDirection ftCase)
+{
+ std::vector<complex_t> cdata;
+ cdata.resize(data.size());
+ for (size_t i = 0; i < data.size(); i++) {
+ cdata[i] = complex_t(data[i], 0);
}
+ return MathFunctions::FastFourierTransform(cdata, ftCase);
+}
+
+//! convolution of two real vectors of equal size
+
+std::vector<complex_t>
+MathFunctions::ConvolveFFT(const std::vector<double> &signal,
+ const std::vector<double> &resfunc)
+{
+ if (signal.size() != resfunc.size()) {
+ throw std::runtime_error("MathFunctions::ConvolveFFT() -> This convolution works only for "
+ "two vectors of equal size. Use Convolve class instead.");
+ }
+ std::vector<complex_t> fft_signal
+ = MathFunctions::FastFourierTransform(signal, MathFunctions::FORWARD_FFT);
+ std::vector<complex_t> fft_resfunc
+ = MathFunctions::FastFourierTransform(resfunc, MathFunctions::FORWARD_FFT);
+
+ std::vector<complex_t> fft_prod;
+ fft_prod.resize(fft_signal.size());
+ for (size_t i = 0; i < fft_signal.size(); i++) {
+ fft_prod[i] = fft_signal[i] * fft_resfunc[i];
+ }
+
+ std::vector<complex_t> result
+ = MathFunctions::FastFourierTransform(fft_prod, MathFunctions::BACKWARD_FFT);
+ return result;
+}
+
+
+// ************************************************************************** //
+// Random number generators
+// ************************************************************************** //
+
+double MathFunctions::GenerateUniformRandom()
+{
+ int random_int = std::rand();
+ return (double)random_int / RAND_MAX;
+}
+
+double MathFunctions::GenerateNormalRandom(double average, double std_dev)
+{
+ return GenerateStandardNormalRandom()*std_dev + average;
+}
+
+double MathFunctions::GenerateStandardNormalRandom() // using GSL
+{
+ gsl_rng *r;
+ r = gsl_rng_alloc(gsl_rng_ranlxs2);
+ double result = gsl_ran_ugaussian(r);
+ gsl_rng_free(r);
return result;
}
diff --git a/Tests/UnitTests/TestCore/MatrixFunctionsTest.h b/Tests/UnitTests/TestCore/MatrixFunctionsTest.h
deleted file mode 100644
index 17c251235d3..00000000000
--- a/Tests/UnitTests/TestCore/MatrixFunctionsTest.h
+++ /dev/null
@@ -1,77 +0,0 @@
-#ifndef MATRIXFUNCTIONSTEST_H
-#define MATRIXFUNCTIONSTEST_H
-
-#include "MathFunctions.h"
-#include "gtest/gtest.h"
-
-class MatrixFunctionsTest : public::testing::Test
-{
-protected:
- MatrixFunctionsTest();
- virtual ~MatrixFunctionsTest(){}
-
- Eigen::Matrix2cd matrix_2cd;
-public:
- EIGEN_MAKE_ALIGNED_OPERATOR_NEW
-};
-
-MatrixFunctionsTest::MatrixFunctionsTest(){
- matrix_2cd(0,0) = complex_t(2.0, 1.0);
- matrix_2cd(0,1) = complex_t(3.0, 1.0);
- matrix_2cd(1,0) = complex_t(4.0, 1.0);
- matrix_2cd(1,1) = complex_t(5.0, 1.0);
-}
-
-TEST_F(MatrixFunctionsTest, Norm)
-{
-
- Eigen::Matrix2d matrix_norm = MathFunctions::Norm(matrix_2cd);
- EXPECT_EQ(std::norm((complex_t)matrix_2cd(0,0)), matrix_norm(0,0));
- EXPECT_EQ(std::norm((complex_t)matrix_2cd(0,1)), matrix_norm(0,1));
- EXPECT_EQ(std::norm((complex_t)matrix_2cd(1,0)), matrix_norm(1,0));
- EXPECT_EQ(std::norm((complex_t)matrix_2cd(1,1)), matrix_norm(1,1));
-}
-
-TEST_F(MatrixFunctionsTest, Abs)
-{
- Eigen::Matrix2d matrix_abs = MathFunctions::Abs(matrix_2cd);
- EXPECT_EQ(std::abs((complex_t)matrix_2cd(0,0)), matrix_abs(0,0));
- EXPECT_EQ(std::abs((complex_t)matrix_2cd(0,1)), matrix_abs(0,1));
- EXPECT_EQ(std::abs((complex_t)matrix_2cd(1,0)), matrix_abs(1,0));
- EXPECT_EQ(std::abs((complex_t)matrix_2cd(1,1)), matrix_abs(1,1));
-}
-
-TEST_F(MatrixFunctionsTest, ProductByElement)
-{
- Eigen::Matrix2cd matrix_2cd2;
- matrix_2cd2(0,0) = complex_t(2.0, 1.0);
- matrix_2cd2(0,1) = complex_t(3.0, 1.0);
- matrix_2cd2(1,0) = complex_t(4.0, 1.0);
- matrix_2cd2(1,1) = complex_t(5.0, 1.0);
-
- Eigen::Matrix2cd matrix_prod = MathFunctions::ProductByElement(matrix_2cd, matrix_2cd2);
- EXPECT_EQ(complex_t(3.0, 4.0), matrix_prod(0,0));
- EXPECT_EQ(complex_t(8.0, 6.0), matrix_prod(0,1));
- EXPECT_EQ(complex_t(15.0, 8.0), matrix_prod(1,0));
- EXPECT_EQ(complex_t(24.0, 10.0), matrix_prod(1,1));
-}
-
-TEST_F(MatrixFunctionsTest, Conj)
-{
- Eigen::Matrix2cd matrix_conj = MathFunctions::Conj(matrix_2cd);
- EXPECT_EQ(std::conj((complex_t)matrix_2cd(0,0)), matrix_conj(0,0));
- EXPECT_EQ(std::conj((complex_t)matrix_2cd(0,1)), matrix_conj(0,1));
- EXPECT_EQ(std::conj((complex_t)matrix_2cd(1,0)), matrix_conj(1,0));
- EXPECT_EQ(std::conj((complex_t)matrix_2cd(1,1)), matrix_conj(1,1));
-}
-
-TEST_F(MatrixFunctionsTest, Real)
-{
- Eigen::Matrix2d matrix_conj = MathFunctions::Real(matrix_2cd);
- EXPECT_EQ(((complex_t)matrix_2cd(0,0)).real(), matrix_conj(0,0));
- EXPECT_EQ(((complex_t)matrix_2cd(0,1)).real(), matrix_conj(0,1));
- EXPECT_EQ(((complex_t)matrix_2cd(1,0)).real(), matrix_conj(1,0));
- EXPECT_EQ(((complex_t)matrix_2cd(1,1)).real(), matrix_conj(1,1));
-}
-
-#endif //MATRIXFUNCTIONSTEST_H
diff --git a/Tests/UnitTests/TestCore/SpecialFunctionsTest.h b/Tests/UnitTests/TestCore/SpecialFunctionsTest.h
index 5df36be699f..178b084c53a 100644
--- a/Tests/UnitTests/TestCore/SpecialFunctionsTest.h
+++ b/Tests/UnitTests/TestCore/SpecialFunctionsTest.h
@@ -66,28 +66,28 @@ TEST_F(SpecialFunctionsTest, csinc)
double ph = 2*M_PI*i/24;
//std::cout << "---------------------------------------------------------------------\n";
//std::cout << "phase = " << ph << "\n";
- EXPECT_EQ( MathFunctions::Sinc(complex_t(0,0)), complex_t(1.,0.) );
+ EXPECT_EQ( MathFunctions::sinc(complex_t(0,0)), complex_t(1.,0.) );
complex_t z;
- z=std::polar(1e-17,ph); EXPECT_CNEAR( MathFunctions::Sinc(z), complex_t(1.,0.), eps );
- z=std::polar(2e-17,ph); EXPECT_CNEAR( MathFunctions::Sinc(z), complex_t(1.,0.), eps );
- z=std::polar(5e-17,ph); EXPECT_CNEAR( MathFunctions::Sinc(z), complex_t(1.,0.), eps );
- z=std::polar(1e-16,ph); EXPECT_CNEAR( MathFunctions::Sinc(z), complex_t(1.,0.), eps );
- z=std::polar(2e-16,ph); EXPECT_CNEAR( MathFunctions::Sinc(z), complex_t(1.,0.), eps );
- z=std::polar(5e-16,ph); EXPECT_CNEAR( MathFunctions::Sinc(z), complex_t(1.,0.), eps );
- z=std::polar(1e-15,ph); EXPECT_CNEAR( MathFunctions::Sinc(z), complex_t(1.,0.), eps );
- z=std::polar(1e-13,ph); EXPECT_CNEAR( MathFunctions::Sinc(z), complex_t(1.,0.), eps );
- z=std::polar(1e-11,ph); EXPECT_CNEAR( MathFunctions::Sinc(z), complex_t(1.,0.), eps );
- z=std::polar(1e-9, ph); EXPECT_CNEAR( MathFunctions::Sinc(z), complex_t(1.,0.), eps );
- z=std::polar(5e-8,ph); EXPECT_CNEAR( MathFunctions::Sinc(z), 1.-z*z/6., eps );
- z=std::polar(2e-8,ph); EXPECT_CNEAR( MathFunctions::Sinc(z), 1.-z*z/6., eps );
- z=std::polar(1e-8,ph); EXPECT_CNEAR( MathFunctions::Sinc(z), 1.-z*z/6., eps );
- z=std::polar(5e-7,ph); EXPECT_CNEAR( MathFunctions::Sinc(z), 1.-z*z/6., eps );
- z=std::polar(2e-7,ph); EXPECT_CNEAR( MathFunctions::Sinc(z), 1.-z*z/6., eps );
- z=std::polar(1e-7,ph); EXPECT_CNEAR( MathFunctions::Sinc(z), 1.-z*z/6., eps );
- z=std::polar(1e-6,ph); EXPECT_CNEAR( MathFunctions::Sinc(z), 1.-z*z/6., eps );
- z=std::polar(1e-5,ph); EXPECT_CNEAR( MathFunctions::Sinc(z), 1.-z*z/6., eps );
- z=std::polar(1e-4,ph); EXPECT_CNEAR( MathFunctions::Sinc(z), 1.-z*z/6.*(1.-z*z/20.), eps );
- z=std::polar(1e-3,ph); EXPECT_CNEAR( MathFunctions::Sinc(z), 1.-z*z/6.*(1.-z*z/20.), eps );
+ z=std::polar(1e-17,ph); EXPECT_CNEAR( MathFunctions::sinc(z), complex_t(1.,0.), eps );
+ z=std::polar(2e-17,ph); EXPECT_CNEAR( MathFunctions::sinc(z), complex_t(1.,0.), eps );
+ z=std::polar(5e-17,ph); EXPECT_CNEAR( MathFunctions::sinc(z), complex_t(1.,0.), eps );
+ z=std::polar(1e-16,ph); EXPECT_CNEAR( MathFunctions::sinc(z), complex_t(1.,0.), eps );
+ z=std::polar(2e-16,ph); EXPECT_CNEAR( MathFunctions::sinc(z), complex_t(1.,0.), eps );
+ z=std::polar(5e-16,ph); EXPECT_CNEAR( MathFunctions::sinc(z), complex_t(1.,0.), eps );
+ z=std::polar(1e-15,ph); EXPECT_CNEAR( MathFunctions::sinc(z), complex_t(1.,0.), eps );
+ z=std::polar(1e-13,ph); EXPECT_CNEAR( MathFunctions::sinc(z), complex_t(1.,0.), eps );
+ z=std::polar(1e-11,ph); EXPECT_CNEAR( MathFunctions::sinc(z), complex_t(1.,0.), eps );
+ z=std::polar(1e-9, ph); EXPECT_CNEAR( MathFunctions::sinc(z), complex_t(1.,0.), eps );
+ z=std::polar(5e-8,ph); EXPECT_CNEAR( MathFunctions::sinc(z), 1.-z*z/6., eps );
+ z=std::polar(2e-8,ph); EXPECT_CNEAR( MathFunctions::sinc(z), 1.-z*z/6., eps );
+ z=std::polar(1e-8,ph); EXPECT_CNEAR( MathFunctions::sinc(z), 1.-z*z/6., eps );
+ z=std::polar(5e-7,ph); EXPECT_CNEAR( MathFunctions::sinc(z), 1.-z*z/6., eps );
+ z=std::polar(2e-7,ph); EXPECT_CNEAR( MathFunctions::sinc(z), 1.-z*z/6., eps );
+ z=std::polar(1e-7,ph); EXPECT_CNEAR( MathFunctions::sinc(z), 1.-z*z/6., eps );
+ z=std::polar(1e-6,ph); EXPECT_CNEAR( MathFunctions::sinc(z), 1.-z*z/6., eps );
+ z=std::polar(1e-5,ph); EXPECT_CNEAR( MathFunctions::sinc(z), 1.-z*z/6., eps );
+ z=std::polar(1e-4,ph); EXPECT_CNEAR( MathFunctions::sinc(z), 1.-z*z/6.*(1.-z*z/20.), eps );
+ z=std::polar(1e-3,ph); EXPECT_CNEAR( MathFunctions::sinc(z), 1.-z*z/6.*(1.-z*z/20.), eps );
}
}
diff --git a/Tests/UnitTests/TestCore/main.cpp b/Tests/UnitTests/TestCore/main.cpp
index 37bc2510515..c39b05f3a41 100644
--- a/Tests/UnitTests/TestCore/main.cpp
+++ b/Tests/UnitTests/TestCore/main.cpp
@@ -39,7 +39,6 @@
#include "MatrixRTCoefficientsTest.h"
#include "ScalarSpecularInfoMapTest.h"
#include "MatrixSpecularInfoMapTest.h"
-#include "MatrixFunctionsTest.h"
#include "FixedBinAxisTest.h"
#include "VariableBinAxisTest.h"
#include "ConstKBinAxisTest.h"
@@ -84,4 +83,3 @@ int main(int argc, char** argv)
// run all google tests
return RUN_ALL_TESTS();
}
-
--
GitLab