Commit de4cabd9 authored by Sebastian Wolf's avatar Sebastian Wolf

Some cleanup in UQ

parent 8394d019
Pipeline #1099 passed with stages
in 15 minutes and 54 seconds
......@@ -29,7 +29,7 @@
namespace xerus {
/**
* @brief Base class for all xerus optimization algorithms, allowing a uniform set of settings
* @brief Base class for (in future) all xerus optimization algorithms, allowing a uniform set of settings.
*/
class OptimizationAlgorithm {
public:
......@@ -42,12 +42,15 @@ namespace xerus {
///@brief The target residual norm at which the algorithm shall stop.
double targetRelativeResidual;
///@brief Minimal relative decrease of the residual norm ( newRes/oldRes ) until either the ranks are increased (if allowed) or the algorithm stops.
///@brief Minimal decrease of the residual norm ( newRes/oldRes ) until either the ranks are increased (if allowed) or the algorithm stops.
double minimalResidualNormDecrease;
///@brief Number of iterations used to check for stopping criteria (e.g. residual[iterations] <= residual[iteration-tracking]*pow(minimalResidualNormDecrease, tracking) )
size_t tracking = 10;
///@brief PerformanceData object used to record the performance of the algorithm.
PerformanceData perfData;
protected:
OptimizationAlgorithm(const size_t _minIterations, const size_t _maxIterations, const double _targetRelativeResidual, const double _minimalResidualNormDecrease);
......
......@@ -56,6 +56,9 @@ namespace xerus {namespace uq {
///@brief Evaluates the solution for the specified parameters
Tensor evaluate(const TTTensor& _x, const std::vector<double>& _parameters, const PolynomBasis _basisType);
///@brief Calculates the mean of all given samples.
Tensor sample_mean(const std::vector<Tensor>& _samples);
class UQMeasurementSet {
......@@ -74,11 +77,16 @@ namespace xerus {namespace uq {
void clear();
};
Tensor sample_mean(const std::vector<Tensor>& _samples);
TTTensor initial_guess(const Tensor& _mean, const UQMeasurementSet& _measurments, const PolynomBasis _polyBasis, const std::vector<size_t>& _dimensions);
/**
* @brief: Calculates an inital guess for the UQ solution using the mean an linear distortions.
* @param _mean An approximation for the mean (expectation value) of the solution, e.g. by calculating the mean of all samples (see sample_mean).
* @param _parameterVectors The parameterVectors for measurements corresponding to linear distortions, i.e. only one parameter should be non-zero (and small) for each measurement.
* @param _solutions The solutions corresponding to the @a _parameterVectors
* @param _polyBasis The PolynomBasis to use for the UQ formulation.
* @param _dimensions The dimensions of the solutions tensor.
*/
TTTensor initial_guess(const Tensor& _mean, const std::vector<std::vector<double>>& _parameterVectors, const std::vector<Tensor>& _solutions, const PolynomBasis _polyBasis, const std::vector<size_t>& _dimensions);
}}
......@@ -29,17 +29,40 @@
namespace xerus { namespace uq {
void uq_adf(TTTensor& _x, const UQMeasurementSet& _measurments, const PolynomBasis _basisType, const double _targetEps = 1e-8, const size_t _maxItr = 0);
///@brief Inplace variant of the UQ ADF to find a solution @a _x (with the inital dimensions and rank) for a given set of @a _measurements.
void uq_adf(TTTensor& _x, const UQMeasurementSet& _measurements, const PolynomBasis _basisType, const double _targetEps = 1e-8, const size_t _maxItr = 0);
TTTensor uq_adf(const UQMeasurementSet& _initalMeasurments, const UQMeasurementSet& _measurments, const PolynomBasis _basisType, const std::vector<size_t>& _dimensions, const double _targetEps = 1e-8, const size_t _maxItr = 0);
///@brief Inplace variant of the rank-adaptive UQ ADF, to find a solution @a _x (with the inital dimensions) for a given set of @a _measurements.
void uq_ra_adf(TTTensor& _x, const UQMeasurementSet& _measurements, const PolynomBasis _basisType, const double _targetEps = 1e-8, const size_t _maxItr = 0, const double _initalRankEps = 1e-2);
void uq_ra_adf(TTTensor& _x, const UQMeasurementSet& _measurments, const PolynomBasis _basisType, const double _targetEps = 1e-8, const size_t _maxItr = 0, const double _initalRankEps = 1e-2);
TTTensor uq_ra_adf(const UQMeasurementSet& _measurments, const PolynomBasis _basisType, const std::vector<size_t>& _dimensions, const double _targetEps = 1e-8, const size_t _maxItr = 0);
///@brief Rank-adaptive UQ ADF to find a solution with given @a _dimensions for a given set of @a _measurements.
TTTensor uq_ra_adf(const UQMeasurementSet& _measurements, const PolynomBasis _basisType, const std::vector<size_t>& _dimensions, const double _targetEps = 1e-8, const size_t _maxItr = 0);
/**
* @brief Rank adaptive ADF to calculate the UQ solution for given measurements.
* @param _positions The positions of the measurements.
* @param _solutions The measured solutions corresponding to the @a _positions.
* @param _dimensions The dimensions of the final solution tensor.
* @param _targetEps TODO describe effect
* @param _maxItr Maximal number of iterations to perform.
*/
TTTensor uq_ra_adf(const std::vector<std::vector<Tensor>>& _positions, const std::vector<Tensor>& _solutions, const std::vector<size_t>& _dimensions, const double _targetEps = 1e-8, const size_t _maxItr = 0);
/**
* @brief Rank adaptive ADF to calculate the UQ solution for given weighted measurements.
* @param _positions The positions of the measurements.
* @param _solutions The measured solutions corresponding to the @a _positions.
* @param _weights Weights for the individual measurements.
* @param _dimensions The dimensions of the final solution tensor.
* @param _targetEps TODO describe effect
* @param _maxItr Maximal number of iterations to perform.
*/
TTTensor uq_ra_adf(const std::vector<std::vector<Tensor>>& _positions, const std::vector<Tensor>& _solutions, const std::vector<double>& _weights, const std::vector<size_t>& _dimensions, const double _targetEps = 1e-8, const size_t _maxItr = 0);
TTTensor uq_ra_adf_iv(TTTensor& _x, const UQMeasurementSet& _measurments, const PolynomBasis _basisType, const double _targetEps = 1e-8, const size_t _maxItr = 0);
// TODO How is this different to the inplace variant? Also very confusing that it changes _x inplace AND returns it.
TTTensor uq_ra_adf_iv(TTTensor& _x, const UQMeasurementSet& _measurements, const PolynomBasis _basisType, const double _targetEps = 1e-8, const size_t _maxItr = 0);
}}
......@@ -34,13 +34,13 @@
namespace xerus {
template<class PositionType> class MeasurementSet {
public:
///@brief Vector containing the positions of the measurments.
///@brief Vector containing the positions of the measurements.
std::vector<std::vector<PositionType>> positions;
///@brief Vector containing the measured values.
std::vector<value_t> measuredValues;
///@brief If empty no weights are considered, otherwise a vector containign the weights for the measurments.
///@brief If empty no weights are considered, otherwise a vector containign the weights for the measurements.
std::vector<value_t> weights;
......@@ -60,8 +60,11 @@ namespace xerus {
///@brief Returns the order of the tensor that is measured.
size_t order() const;
///@brief Returns the 2-norm of the measurments, with weights considered if appropriate.
///@brief Returns the 2-norm of the measurements, with weights considered if appropriate.
value_t norm_2() const;
///@brief Removes all measurements from the set.
void clear();
///@brief Add a measurment at @a _position with @a _value to the set.
void add(const std::vector<PositionType>& _position, const value_t _measuredValue);
......@@ -69,10 +72,10 @@ namespace xerus {
///@brief Add a measurment at @a _position with @a _value and @a _weight to the set.
void add(const std::vector<PositionType>& _position, const value_t _measuredValue, const value_t _weight);
///@brief Sort the measurments.
///@brief Sort the measurements.
void sort();
///@brief Add noise with relative 2-norm @a _epsilon to the measurments.
///@brief Add noise with relative 2-norm @a _epsilon to the measurements.
void add_noise(const double _epsilon);
///@brief Set the measuredValues equal to the ones measured from @a _solution.
......@@ -111,7 +114,7 @@ namespace xerus {
};
/**
* @brief Class used to represent a set of single point measurments.
* @brief Class used to represent a set of single point measurements.
*/
class SinglePointMeasurementSet : public MeasurementSet<size_t> {
public:
......@@ -145,7 +148,7 @@ namespace xerus {
};
/**
* @brief Class used to represent a set of rank-one measurments.
* @brief Class used to represent a set of rank-one measurements.
*/
class RankOneMeasurementSet : public MeasurementSet<Tensor> {
public:
......
......@@ -26,5 +26,11 @@
namespace xerus {
OptimizationAlgorithm::OptimizationAlgorithm(const size_t _minIterations, const size_t _maxIterations, const double _targetRelativeResidual, const double _minimalResidualNormDecrease) : minIterations(_minIterations), maxIterations(_maxIterations), targetRelativeResidual(_targetRelativeResidual), minimalResidualNormDecrease(_minimalResidualNormDecrease) {}
OptimizationAlgorithm::OptimizationAlgorithm(const size_t _minIterations, const size_t _maxIterations, const double _targetRelativeResidual, const double _minimalResidualNormDecrease) :
minIterations(_minIterations),
maxIterations(_maxIterations),
targetRelativeResidual(_targetRelativeResidual),
minimalResidualNormDecrease(_minimalResidualNormDecrease),
perfData(false, false)
{}
} // namespace xerus
......@@ -164,6 +164,19 @@ namespace xerus { namespace uq {
}
Tensor sample_mean(const std::vector<Tensor>& _samples) {
REQUIRE(_samples.size() > 0, "Need at least one measurment.");
// Calc mean
Tensor mean(_samples.front().dimensions);
for(const auto& samp : _samples) {
mean += samp;
}
mean /= double(_samples.size());
return mean;
}
void UQMeasurementSet::add(const std::vector<double>& _paramVec, const Tensor& _solution) {
parameterVectors.push_back(_paramVec);
solutions.push_back(_solution);
......@@ -181,23 +194,11 @@ namespace xerus { namespace uq {
}
Tensor sample_mean(const std::vector<Tensor>& _samples) {
REQUIRE(_samples.size() > 0, "Need at least one measurment.");
// Calc mean
Tensor mean({_samples.front().size});
for(const auto& samp : _samples) {
mean += samp;
}
mean /= double(_samples.size());
return mean;
}
TTTensor initial_guess(const Tensor& _mean, const UQMeasurementSet& _measurments, const PolynomBasis _polyBasis, const std::vector<size_t>& _dimensions) {
REQUIRE(_measurments.parameterVectors.size() > 0, "Need at least one measurment.");
REQUIRE(_measurments.parameterVectors.size() == _measurments.solutions.size(), "Invalid measurments.");
REQUIRE(_dimensions.front() == _measurments.solutions.front().dimensions.front(), "Inconsitend measurments and dimensions.");
// initial_guess(const Tensor& _mean, const UQMeasurementSet& _measurments, const PolynomBasis _polyBasis, const std::vector<size_t>& _dimensions)
TTTensor initial_guess(const Tensor& _mean, const std::vector<std::vector<double>>& _parameterVectors, const std::vector<Tensor>& _solutions, const PolynomBasis _polyBasis, const std::vector<size_t>& _dimensions){
REQUIRE(_parameterVectors.size() > 0, "Need at least one measurment.");
REQUIRE(_parameterVectors.size() == _solutions.size(), "Invalid measurments.");
REQUIRE(_dimensions.front() == _solutions.front().dimensions.front(), "Inconsitend measurments and dimensions.");
TTTensor initalGuess(_dimensions);
Tensor mean = _mean;
......@@ -214,9 +215,9 @@ namespace xerus { namespace uq {
// Calc linear terms
std::set<size_t> usedParams;
for(size_t m = 0; m < _measurments.parameterVectors.size(); ++m) {
const auto& paramVec = _measurments.parameterVectors[m];
const auto& sol = _measurments.solutions[m];
for(size_t m = 0; m < _parameterVectors.size(); ++m) {
const auto& paramVec = _parameterVectors[m];
const auto& sol = _solutions[m];
REQUIRE(paramVec.size()+1 == initalGuess.degree(), "Invalid parameter vector");
......@@ -226,12 +227,12 @@ namespace xerus { namespace uq {
for(size_t i = 0; i < paramVec.size(); ++i) {
if(std::abs(paramVec[i]) > 0.0) {
if(misc::contains(usedParams, i)) {
LOG(info, "Skipping douplicate parameter " << i);
LOG(info, "Skipping duplicate parameter " << i);
skip = true;
continue;
}
REQUIRE(p == initalGuess.degree(), "Parameters contains several non-zero entries: " << paramVec);
REQUIRE(!misc::contains(usedParams, i), "Parameters " << i << " appears twice!" << _measurments.parameterVectors);
REQUIRE(!misc::contains(usedParams, i), "Parameters " << i << " appears twice!" << _parameterVectors);
usedParams.emplace(i);
p = i;
}
......@@ -256,7 +257,7 @@ namespace xerus { namespace uq {
}
initalGuess.round(1e-4);
LOG(UQ_Inital_Guess, "Found linear terms for " << usedParams << ". Post roundign ranks: " << initalGuess.ranks());
LOG(UQ_Inital_Guess, "Found linear terms for " << usedParams << ". Post rounding ranks: " << initalGuess.ranks());
return initalGuess;
}
......
......@@ -541,18 +541,6 @@ namespace xerus { namespace uq { namespace impl_uqRaAdf {
}
TTTensor uq_adf(const UQMeasurementSet& _initalMeasurments, const UQMeasurementSet& _measurments, const PolynomBasis _basisType, const std::vector<size_t>& _dimensions, const double _targetEps, const size_t _maxItr) {
REQUIRE(_measurments.parameterVectors.size() == _measurments.solutions.size(), "Invalid measurments");
REQUIRE(_dimensions.front() == _measurments.solutions.front().size, "Inconsitent spacial dimension");
TTTensor x = initial_guess(sample_mean(_measurments.solutions), _initalMeasurments, _basisType, _dimensions);
impl_uqRaAdf::InternalSolver<1> solver(x, _measurments, _basisType, _maxItr, _targetEps, 0.0);
solver.solve();
return x;
}
void uq_ra_adf(TTTensor& _x, const UQMeasurementSet& _measurments, const PolynomBasis _basisType, const double _targetEps, const size_t _maxItr, const double _initalRankEps) {
REQUIRE(_measurments.parameterVectors.size() == _measurments.solutions.size(), "Invalid measurments");
REQUIRE(_x.dimensions.front() == _measurments.solutions.front().size, "Inconsitent spacial dimension");
......@@ -584,6 +572,7 @@ namespace xerus { namespace uq { namespace impl_uqRaAdf {
solver.solve();
return x;
}
TTTensor uq_ra_adf(const std::vector<std::vector<Tensor>>& _positions, const std::vector<Tensor>& _solutions, const std::vector<size_t>& _dimensions, const double _targetEps, const size_t _maxItr) {
REQUIRE(_positions.size() == _solutions.size(), "Invalid measurments");
......@@ -607,6 +596,7 @@ namespace xerus { namespace uq { namespace impl_uqRaAdf {
solver.solve();
return x;
}
TTTensor uq_ra_adf(const std::vector<std::vector<Tensor>>& _positions, const std::vector<Tensor>& _solutions, const std::vector<double>& _weights, const std::vector<size_t>& _dimensions, const double _targetEps, const size_t _maxItr) {
REQUIRE(_positions.size() == _solutions.size(), "Invalid measurments");
......@@ -630,6 +620,7 @@ namespace xerus { namespace uq { namespace impl_uqRaAdf {
solver.solve();
return x;
}
TTTensor uq_ra_adf_iv(TTTensor& _x, const UQMeasurementSet& _measurments, const PolynomBasis _basisType, const double _targetEps, const size_t _maxItr) {
REQUIRE(_measurments.parameterVectors.size() == _measurments.solutions.size(), "Invalid measurments");
......
......@@ -94,6 +94,14 @@ namespace xerus {
}
template<class PositionType>
void MeasurementSet<PositionType>::clear() {
positions.clear();
measuredValues.clear();
weights.clear();
}
template<class PositionType>
void MeasurementSet<PositionType>::add(const std::vector<PositionType>& _position, const value_t _measuredValue) {
XERUS_REQUIRE_TEST;
......
Markdown is supported
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