Commit 0b186260 authored by Nabil Chouika's avatar Nabil Chouika
Browse files

refs#16

In ElementaryUtils:
- Use command "executable" in the Ant scripts instead of the deprecated command "command".

In NumA++:
- Overload functors to allow call without parameters.
- Restructuring of folders for Functors.
- Rename LinearSystem to LinAlgUtils and add rank routine.
- New EigenUtils to encapsulate conversions to or from Eigen. In order for PARTONS to be independent of Eigen.
- New methods in MatrixD (appendLine, etc).
- New constructors for VectorD.
- New static method in Interval to compute nodes.

In NumA++_Test:
- New tests for Eigen.

In PARTONS:
- New modules IncompleteGPDModule and RadonInverseModule.
- First method for Radon inversion: Random brute force with constant piecewise functions (RandBFConstPW). Still no inversion, only matrix construction.
- Addition of DDGauge type and modification of DDType.
- New OverlapMMR15 model.
parent a6c21879
......@@ -30,6 +30,7 @@ public:
virtual double operator()(VectorD &variables,
std::vector<double> &parameters) = 0;
virtual double operator()(VectorD &variables) = 0;
};
} // namespace NumA
......
......@@ -35,6 +35,10 @@ public:
std::vector<double> &parameters) {
return ((*m_pObject).*m_pFunction)(variables, parameters);
}
inline double operator()(VectorD &variables) {
std::vector<double> parameters;
return ((*m_pObject).*m_pFunction)(variables, parameters);
}
private:
PointerToObj* m_pObject;
......
......@@ -27,6 +27,7 @@ public:
}
virtual double operator()(double x, std::vector<double> &parameters) = 0;
virtual double operator()(double x) = 0;
};
} // namespace NumA
......
......@@ -33,6 +33,10 @@ public:
inline double operator()(double x, std::vector<double> &parameters) {
return ((*m_pObject).*m_pFunction)(x, parameters);
}
inline double operator()(double x) {
std::vector<double> parameters;
return ((*m_pObject).*m_pFunction)(x, parameters);
}
private:
PointerToObj* m_pObject;
......
......@@ -11,9 +11,9 @@
#include <string>
#include <vector>
#include "../../functor/one_dimension/Functor1D.h"
#include "../../utils/Errors.h"
#include "../../utils/Tolerances.h"
#include "Functor1D.h"
#include "IntegratorType1D.h"
namespace ElemUtils {
......
......@@ -8,16 +8,16 @@
* @version 1.0
*/
#include <stddef.h>
#include <stdlib.h>
#include <cmath>
#include <iostream>
#include <limits>
#include <vector>
#include <stddef.h>
#include <stdlib.h>
#include "../../../utils/Errors.h"
#include "../../../utils/Tolerances.h"
#include "../FunctionType1D.h"
#include "../../../functor/one_dimension/FunctionType1D.h"
//#include "../../../utils/Errors.h"
//#include "../../../utils/Tolerances.h"
#include "../Integrator1D.h"
#define GSL_DBL_EPSILON 2.2204460492503131e-16
......
/**
* @file EigenUtils.h
* @author Nabil Chouika (Irfu/SPhN, CEA Saclay)
* @date 11 oct. 2016
* @version 1.0
*/
#ifndef EIGENUTILS_H_
#define EIGENUTILS_H_
#include <Eigen/Core>
#include "../matrix/MatrixD.h"
#include "../vector/VectorD.h"
namespace NumA {
/**
* @class EigenUtils
* @brief Tools for the Eigen wrapper.
*/
class EigenUtils {
public:
static Eigen::MatrixXd convertToEigen(const MatrixD & A);
static Eigen::VectorXd convertToEigen(const VectorD & V);
static VectorD convertToVectorD(const Eigen::VectorXd & V);
};
} /* namespace NumA */
#endif /* EIGENUTILS_H_ */
/**
* @file LinearSystem.h
* @file LinAlgUtils.h
* @author Nabil CHOUIKA (SPhN / CEA Saclay)
* @date Feb 4 2016
* @version 1.0
*/
#ifndef LINEARSYSTEM_H_
#define LINEARSYSTEM_H_
#ifndef LINALGUTILS_H_
#define LINALGUTILS_H_
#include "../matrix/MatrixD.h"
#include "../vector/VectorD.h"
namespace NumA {
......@@ -17,11 +18,11 @@ class MatrixD;
namespace NumA {
/**
* @class LinearSystem
* @brief Linear system solver.
* @class LinAlgUtils
* @brief Linear algebra routines such as linear solvers. Serves as wrapper for Eigen.
*/
class LinearSystem {
class LinAlgUtils {
public:
enum Method {
ColQR = 0, FullQR, QR, PartialLU, FullLU
......@@ -35,8 +36,16 @@ public:
* @return X Solution: VectorD
*/
static VectorD solve(const MatrixD & A, const VectorD & Y, Method method = ColQR);
/**
* @brief Computes the rank of the matrix with the given method.
* @param A MatrixD
* @param method
* @return rank unsigned int
*/
static unsigned int rank(const MatrixD & A, Method method = ColQR);
};
} // namespace NumA
#endif /* LINEARSYSTEM_H_ */
#endif /* LINALGUTILS_H_ */
......@@ -83,7 +83,9 @@ public:
void assign(const size_t _rowsNumber, const size_t _columnsNumber,
double value = 0.);
void setLine(unsigned int i, const NumA::VectorD& line);
void setLine(size_t i, const NumA::VectorD& line);
void addLine(size_t i, const NumA::VectorD& line);
void appendLine(const NumA::VectorD& line);
NumA::VectorD getLine(const size_t lineIndex) const;
// Matrix/vector operations
......
......@@ -11,7 +11,14 @@
#include <string>
#include <vector>
#include <stddef.h>
#include "../matrix/MatrixD.h"
//#include "../matrix/MatrixD.h"
namespace NumA {
class Vector2D;
class Vector3D;
class Vector4D;
} /* namespace NumA */
namespace NumA {
......@@ -49,6 +56,30 @@ public:
*/
VectorD(const std::vector<double> &vector);
/**
* Copy constructor.
* @param vector VectorD
*/
VectorD(const VectorD &vector);
/**
* Copy constructor.
* @param vector Vector2D
*/
VectorD(const Vector2D &vector);
/**
* Copy constructor.
* @param vector Vector3D
*/
VectorD(const Vector3D &vector);
/**
* Copy constructor.
* @param vector Vector4D
*/
VectorD(const Vector4D &vector);
/**
* Default destructor
*/
......
......@@ -8,11 +8,11 @@
* @version 1.0
*/
#include <stddef.h>
#include <cmath>
#include <iostream>
#include <stdexcept>
#include <vector>
#include <stddef.h>
namespace NumA {
......@@ -64,7 +64,7 @@ public:
// Newton step
x = x0 - result / result_deriv;
noConvergence = std::fabs(x - x0) > absTolerance;
noConvergence = fabs(x - x0) > absTolerance;
x0 = x;
}
......@@ -97,7 +97,7 @@ public:
x = x1 - v1 * (x1 - x0) / (v1 - v0);
}
noConvergence = std::fabs(x - x1) > absTolerance;
noConvergence = fabs(x - x1) > absTolerance;
x0 = x1;
v0 = v1;
......
......@@ -10,7 +10,7 @@
#include <stddef.h>
#include "../linear_algebra/linear_system/LinearSystem.h"
#include "../linear_algebra/eigen/LinAlgUtils.h"
#include "../linear_algebra/matrix/MatrixD.h"
#include "../linear_algebra/vector/VectorD.h"
......@@ -32,13 +32,13 @@ namespace NumA {
*/
class NewtonMD {
public:
NewtonMD(int N = 1, LinearSystem::Method linMethod = LinearSystem::ColQR);
NewtonMD(int N = 1, LinAlgUtils::Method linMethod = LinAlgUtils::ColQR);
virtual ~NewtonMD();
int getN() const;
void setN(int n);
LinearSystem::Method getLinMethod() const;
void setLinMethod(LinearSystem::Method linMethod);
LinAlgUtils::Method getLinMethod() const;
void setLinMethod(LinAlgUtils::Method linMethod);
template<class OBJ, typename FUNC, typename JACOB>
VectorD solve(OBJ* object, FUNC G,
......@@ -81,7 +81,7 @@ public:
private:
int m_N; ///< Dimension of the problem
LinearSystem::Method m_linMethod; ///< Method used for solving linear systems
LinAlgUtils::Method m_linMethod; ///< Method used for solving linear systems
};
} /* namespace NumA */
......
......@@ -8,8 +8,8 @@
* @version 1.0
*/
#include "../integration/many_dimension/FunctorMD.h"
#include "../integration/one_dimension/Functor1D.h"
#include "../functor/multi_dimension/FunctorMD.h"
#include "../functor/one_dimension/Functor1D.h"
namespace NumA {
......@@ -21,7 +21,7 @@ namespace NumA {
class FunctorUtils {
public:
template<typename PointerToObj, typename PointerToMemFn>
static Functor1D<PointerToObj, PointerToMemFn>* newIntegrationFunctor(
static Functor1D<PointerToObj, PointerToMemFn>* newFunctor1D(
PointerToObj* object, PointerToMemFn function) {
return new Functor1D<PointerToObj, PointerToMemFn>(object, function);
}
......
......@@ -8,8 +8,10 @@
* @version 1.0
*/
#include <cmath>
#include <stdexcept>
#include <vector>
#include <cstddef>
namespace NumA {
......@@ -59,6 +61,34 @@ public:
return steps;
}
static std::vector<T> computeNodes(T lowerBound, T upperBound, size_t num =
10, Interval::StepMode stepMode = NORMAL) {
T node = lowerBound;
T step;
std::vector<T> nodes(num, node);
switch (stepMode) {
case NORMAL:
step = (upperBound - lowerBound) / (num - 1);
for (size_t i = 1; i < num; i++) {
node += step;
nodes[i] = node;
}
break;
case LOG:
step = pow((upperBound / lowerBound), 1 / (num - 1));
for (size_t i = 1; i < num; i++) {
node *= step;
nodes[i] = node;
}
break;
default:
throw std::runtime_error(
"(Interval::computeNodes) Unknown step mode");
break;
}
return nodes;
}
T getLowerBound() const {
return m_lowerBound;
}
......
#include "../../../../include/NumA/integration/one_dimension/DExpIntegrator1D.h"
#include <ElementaryUtils/logger/LoggerManager.h>
#include <cmath>
#include <limits>
#include <vector>
#include <ElementaryUtils/logger/LoggerManager.h>
#include "../../../../include/NumA/integration/one_dimension/FunctionType1D.h"
#include "../../../../include/NumA/utils/Errors.h"
#include "../../../../include/NumA/utils/Tolerances.h"
#include "../../../../include/NumA/functor/one_dimension/FunctionType1D.h"
//#include "../../../../include/NumA/utils/Errors.h"
//#include "../../../../include/NumA/utils/Tolerances.h"
namespace NumA {
......
#include "../../../../include/NumA/integration/one_dimension/QuadratureIntegrator1D.h"
#include <string>
#include <vector>
#include <ElementaryUtils/parameters/GenericType.h>
#include <ElementaryUtils/parameters/Parameters.h>
#include "../../../../include/NumA/integration/one_dimension/FunctionType1D.h"
#include "../../../../include/NumA/functor/one_dimension/FunctionType1D.h"
namespace NumA {
......
#include "../../../../include/NumA/linear_algebra/eigen/EigenUtils.h"
namespace NumA {
Eigen::MatrixXd EigenUtils::convertToEigen(const MatrixD& A) {
Eigen::MatrixXd EigenA(A.rows(), A.cols());
for (unsigned int i = 0; i < A.rows(); i++) {
for (unsigned int j = 0; j < A.cols(); j++) {
EigenA(i, j) = A.at(i, j);
}
}
return EigenA;
}
Eigen::VectorXd EigenUtils::convertToEigen(const VectorD& V) {
Eigen::VectorXd EigenV(V.size());
for (unsigned int i = 0; i < V.size(); i++) {
EigenV(i) = V.at(i);
}
return EigenV;
}
VectorD EigenUtils::convertToVectorD(const Eigen::VectorXd& EigenV) {
VectorD V(EigenV.size());
for (unsigned int i = 0; i < V.size(); i++) {
V.at(i) = EigenV(i);
}
return V;
}
} /* namespace NumA */
#include "../../../../include/NumA/linear_algebra/linear_system/LinearSystem.h"
#include "../../../../include/NumA/linear_algebra/eigen/LinAlgUtils.h"
#include <Eigen/Core>
#include <Eigen/LU>
......@@ -6,11 +6,13 @@
#include <sstream>
#include <stdexcept>
#include "../../../../include/NumA/linear_algebra/matrix/MatrixD.h"
#include "../../../../include/NumA/linear_algebra/eigen/EigenUtils.h"
//#include "../../../../include/NumA/linear_algebra/matrix/MatrixD.h"
namespace NumA {
VectorD LinearSystem::solve(const MatrixD & A, const VectorD & Y, Method method) {
VectorD LinAlgUtils::solve(const MatrixD & A, const VectorD & Y, Method method) {
if (A.cols() != A.rows()) {
std::stringstream formatter;
formatter << "[Solve::solve] Matrix column size = " << A.cols()
......@@ -25,20 +27,12 @@ VectorD LinearSystem::solve(const MatrixD & A, const VectorD & Y, Method method)
}
// Converting matrix and vector to Eigen
Eigen::MatrixXd EigenA(A.rows(),A.cols());
for (unsigned int i = 0; i < A.rows(); i++) {
for (unsigned int j = 0; j < A.cols(); j++) {
EigenA(i,j) = A.at(i,j);
}
}
Eigen::VectorXd EigenY(Y.size());
for (unsigned int i = 0; i < Y.size(); i++) {
EigenY(i) = Y.at(i);
}
Eigen::MatrixXd EigenA = EigenUtils::convertToEigen(A);
Eigen::VectorXd EigenY = EigenUtils::convertToEigen(Y);
// Solving the system
Eigen::VectorXd EigenX;
switch(method) {
switch (method) {
case FullQR:
EigenX = EigenA.fullPivHouseholderQr().solve(EigenY);
break;
......@@ -57,12 +51,29 @@ VectorD LinearSystem::solve(const MatrixD & A, const VectorD & Y, Method method)
}
// Converting to VectorD
VectorD X(Y.size());
for (unsigned int i = 0; i < X.size(); i++) {
X.at(i) = EigenX(i);
}
VectorD X = EigenUtils::convertToVectorD(EigenX);
return X;
}
unsigned int LinAlgUtils::rank(const MatrixD& A, Method method) {
// Converting matrix and vector to Eigen
Eigen::MatrixXd EigenA = EigenUtils::convertToEigen(A);
unsigned int rank;
switch (method) {
case FullQR:
rank = EigenA.fullPivHouseholderQr().rank();
break;
case FullLU:
rank = EigenA.fullPivLu().rank();
break;
default:
rank = EigenA.colPivHouseholderQr().rank();
break;
}
return rank;
}
} /* namespace NumA */
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment