Commit f56ae533 authored by Sebastian Wolf's avatar Sebastian Wolf

Work on ASD (use unified interface) and unit tests for it

parent 4b2d5559
......@@ -277,13 +277,13 @@ build/.testObjects/%.o: %.cpp $(MINIMAL_DEPS)
# Build rule for unit test objects
ifdef USE_GCC
build/.unitTestObjects/%.o: %.cxx $(MINIMAL_DEPS) build/.preCompileHeaders/xerus.h.gch
build/.unitTestObjects/%.o: %.cxx $(MINIMAL_DEPS) build/.preCompileHeaders/xerus.h.gch build/.preCompileHeaders/common.hxx.gch
mkdir -p $(dir $@)
$(CXX) -D XERUS_UNITTEST -I build/.preCompileHeaders $< -c $(FLAGS) -MMD -o $@
else
build/.unitTestObjects/%.o: %.cxx $(MINIMAL_DEPS)
mkdir -p $(dir $@)
$(CXX) -D XERUS_UNITTEST -I include $< -c $(FLAGS) -MMD -o $@
$(CXX) -D XERUS_UNITTEST -I include -I src/unitTests $< -c $(FLAGS) -MMD -o $@
endif
......@@ -299,6 +299,12 @@ build/.preCompileHeaders/xerus.h.gch: include/xerus.h $(MINIMAL_DEPS) .git/ORIG_
$(CXX) -D XERUS_UNITTEST $< -c $(FLAGS) -MMD -o $@
# Build rule for the other preCompileHeader
build/.preCompileHeaders/common.hxx.gch: src/unitTests/common.hxx $(MINIMAL_DEPS) .git/ORIG_HEAD
mkdir -p $(dir $@)
$(CXX) -D XERUS_UNITTEST $< -c $(FLAGS) -MMD -o $@
# dummy rule in case files were downloaded without git
.git/ORIG_HEAD:
mkdir -p .git
......
......@@ -47,7 +47,7 @@ namespace xerus {
* @returns the residual @f$|P_\Omega(x-b)|_2@f$ of the final @a _x.
*/
template<class MeasurmentSet>
double operator()(TTTensor& _x, const MeasurmentSet& _measurments, PerformanceData& _perfData) const;
double operator()(TTTensor& _x, const MeasurmentSet& _measurments, PerformanceData& _perfData = NoPerfData) const;
/**
......@@ -59,7 +59,7 @@ namespace xerus {
* @returns the residual @f$|P_\Omega(x-b)|_2@f$ of the final @a _x.
*/
template<class MeasurmentSet>
double operator()(TTTensor& _x, const MeasurmentSet& _measurments, const std::vector<size_t>& _maxRanks, PerformanceData& _perfData) const;
double operator()(TTTensor& _x, const MeasurmentSet& _measurments, const std::vector<size_t>& _maxRanks, PerformanceData& _perfData = NoPerfData) const;
};
/// @brief Default variant of the ADF algorithm
......
......@@ -33,11 +33,13 @@ namespace xerus {
public:
double minRankEps = 1e-4;
double maxRankEps = 0.32;
double epsDecay = 1.1;
double controlSetFraction = 0.1;
double initialRankEps = 5e-3;
double initialRankEps = 4e-3;
/// Basic constructor
ASDVariant(const size_t _maxIterations, const double _targetRelativeResidual, const double _minimalResidualNormDecrease)
......@@ -52,7 +54,7 @@ namespace xerus {
* @param _perfData optinal performanceData object to be used.
* @returns nothing
*/
void operator()(TTTensor& _x, const RankOneMeasurementSet& _measurments, PerformanceData& _perfData) const;
void operator()(TTTensor& _x, const RankOneMeasurementSet& _measurments, PerformanceData& _perfData = NoPerfData) const;
/**
* @brief Tries to reconstruct the (low rank) tensor @a _x from the given measurments.
......@@ -62,7 +64,7 @@ namespace xerus {
* @param _perfData optinal performanceData object to be used.
* @returns nothing
*/
void operator()(TTTensor& _x, const RankOneMeasurementSet& _measurments, const std::vector<size_t>& _maxRanks, PerformanceData& _perfData) const;
void operator()(TTTensor& _x, const RankOneMeasurementSet& _measurments, const std::vector<size_t>& _maxRanks, PerformanceData& _perfData = NoPerfData) const;
};
/// @brief Default variant of the ASD algorithm
......
......@@ -36,5 +36,6 @@ namespace xerus {
class RankOneMeasurementSet;
class PerformanceData;
extern PerformanceData NoPerfData;
} // End xerus namespace
......@@ -34,5 +34,6 @@ XERUS_SET_LOGGING(unit_test, xerus::misc::internal::LOGGING_ON_ERROR)
XERUS_SET_LOGGING(unit_tests, xerus::misc::internal::LOGGING_ON_ERROR)
XERUS_SET_LOGGING(pydebug, xerus::misc::internal::LOGGING_ON_ERROR)
XERUS_SET_LOGGING(ADF, xerus::misc::internal::LOGGING_ON_ERROR)
XERUS_SET_LOGGING(ASD, xerus::misc::internal::LOGGING_ON_ERROR)
XERUS_SET_LOGGING(UQ, xerus::misc::internal::LOGGING_ON_ERROR)
/* */
// Xerus - A General Purpose Tensor Library
// Copyright (C) 2014-2019 Benjamin Huber and Sebastian Wolf.
//
// Xerus is free software: you can redistribute it and/or modify
// it under the terms of the GNU Affero General Public License as published
// by the Free Software Foundation, either version 3 of the License,
// or (at your option) any later version.
//
// Xerus 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 Affero General Public License for more details.
//
// You should have received a copy of the GNU Affero General Public License
// along with Xerus. If not, see <http://www.gnu.org/licenses/>.
//
// For further information on Xerus visit https://libXerus.org
// or contact us at contact@libXerus.org.
#include <common.hxx>
static misc::UnitTest trasd_frr("TRASD", "Fixed_Rank_Recovery", [](){
const std::vector<size_t> orders{2, 3, 5};
const size_t runs = 6;
const size_t successThreshold = 4;
auto myTRASD = TRASD;
myTRASD.maxIterations = 500;
myTRASD.targetRelativeResidual = 1e-6;
for(const auto d : orders) {
size_t residualSucc = 0, ErrorSucc = 0;
for(size_t i = 0; i < runs; ++i) {
const auto dimensions = random_dimensions(d, 1, 8);
if(misc::product(dimensions) == 1) {continue;} // TODO Check why this doesn't work.
const auto ranks = random_low_tt_ranks(dimensions, 1, 6);
const auto target = TTTensor::random(dimensions, ranks);
const double targetNorm = frob_norm(target);
const size_t targetDofs = target.datasize();
// This should allow succesfull reconstruction in almost all cases.
const size_t N = 2*targetDofs+misc::product(target.dimensions)/25;
const auto measurments = RankOneMeasurementSet::random(N, target);
TTTensor solution = TTTensor::random(target.dimensions, ranks);
myTRASD(solution, measurments);
const double error = frob_norm(target-solution)/targetNorm;
const double residual = measurments.test(solution);
if(residual < 1e-6) { residualSucc++; }
if(error < 1e-4) { ErrorSucc++; }
}
MTEST(residualSucc >= successThreshold, "Only " << residualSucc << " of " << runs << " were succesfull in terms of residual.");
MTEST(ErrorSucc >= successThreshold, "Only " << ErrorSucc << " of " << runs << " were succesfull in terms of error.");
}
});
static misc::UnitTest trasd_rar("TRASD", "Rank Adaptive_Recovery", [](){
const std::vector<size_t> orders{2, 3, 4};
const size_t runs = 6;
const size_t successThreshold = 4;
auto myTRASD = TRASD;
myTRASD.maxIterations = 500;
myTRASD.targetRelativeResidual = 1e-6;
for(const auto d : orders) {
size_t residualSucc = 0, ErrorSucc = 0;
for(size_t i = 0; i < runs; ++i) {
const auto dimensions = random_dimensions(d, 1, 7);
if(misc::product(dimensions) == 1) {continue;} // TODO Check why this doesn't work.
const auto ranks = random_low_tt_ranks(dimensions, 1, 5);
const auto target = TTTensor::random(dimensions, ranks);
const double targetNorm = frob_norm(target);
const size_t targetDofs = target.datasize();
// This should allow succesfull reconstruction in almost all cases.
const size_t N = 4*targetDofs+misc::product(target.dimensions)/4;
const auto measurments = RankOneMeasurementSet::random(N, target);
TTTensor solution = TTTensor::random(target.dimensions, std::vector<size_t>(d-1, 1));
myTRASD(solution, measurments, std::vector<size_t>(d-1, 10));
const double error = frob_norm(target-solution)/targetNorm;
const double residual = measurments.test(solution);
if(residual < 1e-4) { residualSucc++; }
if(error < 1e-4) { ErrorSucc++; }
}
MTEST(residualSucc >= successThreshold, "Only " << residualSucc << " of " << runs << " were succesfull in terms of residual.");
MTEST(ErrorSucc >= successThreshold, "Only " << ErrorSucc << " of " << runs << " were succesfull in terms of error.");
}
});
// static misc::UnitTest trasd_frc("TRASD", "Fixed_Rank_Completion", [](){
// const std::vector<size_t> orders{2, 3, 5};
// const size_t runs = 100;
// const size_t successThreshold = 100;
//
// auto myTRASD = TRASD;
//
// myTRASD.maxIterations = 500;
// myTRASD.targetRelativeResidual = 1e-6;
//
// for(const auto d : orders) {
// size_t residualSucc = 0, ErrorSucc = 0;
//
// for(size_t i = 0; i < runs; ++i) {
// const auto dimensions = random_dimensions(d, 1, 8);
//
// if(misc::product(dimensions) == 1) {continue;} // TODO Check why this doesn't work.
//
// const auto ranks = random_low_tt_ranks(dimensions, 1, 6);
//
// const auto target = TTTensor::random(dimensions, ranks);
// const double targetNorm = frob_norm(target);
// const size_t targetDofs = target.datasize();
//
// // This should allow succesfull reconstruction in almost all cases.
// const size_t N = std::min(misc::product(dimensions), 5*targetDofs+misc::product(target.dimensions)/25);
//
// const auto completionMeasurments = RankOneMeasurementSet(SinglePointMeasurementSet::random(N, target), target.dimensions);
//
//
// TTTensor solution = TTTensor::random(target.dimensions, ranks);
// myTRASD(solution, completionMeasurments);
// const double error = frob_norm(target-solution)/targetNorm;
// const double residual = completionMeasurments.test(solution);
// if(residual < 1e-6) { residualSucc++; }
// if(error < 1e-4) { ErrorSucc++; }
// }
//
// MTEST(residualSucc >= successThreshold, "Only " << residualSucc << " of " << runs << " were succesfull in terms of residual.");
// MTEST(ErrorSucc >= successThreshold, "Only " << ErrorSucc << " of " << runs << " were succesfull in terms of error.");
// }
// });
// Xerus - A General Purpose Tensor Library
// Copyright (C) 2014-2019 Benjamin Huber and Sebastian Wolf.
//
// Xerus is free software: you can redistribute it and/or modify
// it under the terms of the GNU Affero General Public License as published
// by the Free Software Foundation, either version 3 of the License,
// or (at your option) any later version.
//
// Xerus 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 Affero General Public License for more details.
//
// You should have received a copy of the GNU Affero General Public License
// along with Xerus. If not, see <http://www.gnu.org/licenses/>.
//
// For further information on Xerus visit https://libXerus.org
// or contact us at contact@libXerus.org.
/**
* @file
* @brief Common functions for the Unittests.
*/
#include "../../include/xerus.h"
#include "../../include/xerus/test/test.h"
#include "../../include/xerus/misc/internal.h"
using namespace xerus;
size_t tt_dofs(const std::vector<size_t>& _dimensions, const std::vector<size_t>& _ranks) {
REQUIRE(_dimensions.size() == _ranks.size()+1, "Inconsitend order.");
const size_t order = _dimensions.size();
size_t dofs = 0;
dofs += _dimensions[0]*_ranks[0];
for(size_t mu = 1; mu+1 < order; ++mu) {
dofs += _ranks[mu-1]*_dimensions[mu]*_ranks[mu];
}
dofs += _ranks[order-1]*_dimensions[order-1];
for(size_t mu = 0; mu+2 < order; ++mu) {
dofs -= misc::sqr(_ranks[mu]);
}
return dofs;
}
std::vector<size_t> random_dimensions(const size_t _order, const size_t _minN = 1, const size_t _maxN = 10) {
std::uniform_int_distribution<size_t> dimDist(_minN, _maxN);
std::vector<size_t> dimensions;
for(size_t d = 0; d < _order; ++d) {
dimensions.push_back(dimDist(misc::randomEngine));
}
return dimensions;
}
std::vector<size_t> random_tt_ranks(const std::vector<size_t>& _dimensions, const size_t _minN = 1, const size_t _maxN = 10) {
std::uniform_int_distribution<size_t> dimDist(_minN, _maxN);
std::vector<size_t> ranks;
for(size_t d = 1; d < _dimensions.size(); ++d) {
ranks.push_back(dimDist(misc::randomEngine));
}
return TTNetwork<false>::reduce_to_maximal_ranks(ranks, _dimensions);
}
std::vector<size_t> random_low_tt_ranks(const std::vector<size_t>& _dimensions, const size_t _minN = 1, size_t _maxN = 10) {
std::vector<size_t> ranks;
do {
std::uniform_int_distribution<size_t> dimDist(_minN, _maxN);
ranks.clear();
for(size_t d = 1; d < _dimensions.size(); ++d) {
ranks.push_back(dimDist(misc::randomEngine));
}
ranks = TTNetwork<false>::reduce_to_maximal_ranks(ranks, _dimensions);
--_maxN;
} while (_maxN > 2 && tt_dofs(_dimensions, ranks) > misc::product(_dimensions)/10);
return ranks;
}
......@@ -42,7 +42,7 @@ namespace xerus { namespace impl_TrASD {
template <size_t P>
class InternalSolver {
class InternalSolver : internal::OptimizationSolver {
private:
///@brief Reference to the external solution (external ownership).
TTTensor& outX;
......@@ -59,35 +59,18 @@ namespace xerus { namespace impl_TrASD {
///@brief Number of measurments (i.e. measurments.size())
const size_t numMeasurments;
///@brief Minimal number of iterations (one iteration = one sweep)
const size_t minIterations;
///@brief Maximal allowed number of iterations (one iteration = one sweep)
const size_t maxIterations;
///@brief The target residual norm at which the algorithm shall stop.
const double targetRelativeResidual;
///@brief Minimal relative decrease of the residual norm ( newRes/oldRes ) until either the ranks are increased (if allowed) or the algorithm stops.
const double minimalResidualNormDecrease;
const size_t tracking;
///@brief Maximally allowed ranks.
const std::vector<size_t> maxRanks;
const double minRankEps;
const double maxRankEps;
const double epsDecay;
const double controlSetFraction;
///@brief: Reference to the performanceData object (external ownership)
PerformanceData& perfData;
///@brief Vector of measurment IDs for each set (0 to P-1). The set P is the control set.
std::vector<std::vector<size_t>> sets;
......@@ -97,7 +80,6 @@ namespace xerus { namespace impl_TrASD {
///@brief L2 Norm of the optimizing sets (0-P-1).
double optNorm;
///@brief Current rankEps
double rankEps;
......@@ -109,8 +91,6 @@ namespace xerus { namespace impl_TrASD {
boost::circular_buffer<std::vector<size_t>> prevRanks;
boost::circular_buffer<double> residuals;
///@brief Stack of pre calculated components from corePosition 0 to d-1
std::vector<std::vector<Tensor>> leftStack;
......@@ -144,12 +124,13 @@ namespace xerus { namespace impl_TrASD {
}
}
InternalSolver( TTTensor& _x,
InternalSolver( const ASDVariant& _optiAlgorithm,
TTTensor& _x,
const RankOneMeasurementSet& _measurments,
const ASDVariant& _optiSettings,
const double _initalRankEps,
const std::vector<size_t>& _maxRanks,
PerformanceData& _perfData ) :
OptimizationSolver(_optiAlgorithm, _perfData),
outX(_x),
x(_x, 0, P),
degree(_x.degree()),
......@@ -157,21 +138,13 @@ namespace xerus { namespace impl_TrASD {
measurments(_measurments),
numMeasurments(_measurments.size()),
minIterations(_optiSettings.minIterations),
maxIterations(_optiSettings.maxIterations),
targetRelativeResidual(_optiSettings.targetRelativeResidual),
minimalResidualNormDecrease(_optiSettings.minimalResidualDecrease),
tracking(_optiSettings.tracking),
maxRanks(TTTensor::reduce_to_maximal_ranks(_maxRanks, _x.dimensions)),
minRankEps(_optiSettings.minRankEps),
epsDecay(_optiSettings.epsDecay),
controlSetFraction(P == 1 ? 0.0 : _optiSettings.controlSetFraction),
perfData(_perfData),
minRankEps(_optiAlgorithm.minRankEps),
maxRankEps(_optiAlgorithm.maxRankEps),
epsDecay(_optiAlgorithm.epsDecay),
controlSetFraction(P == 1 ? 0.0 : _optiAlgorithm.controlSetFraction),
sets(P+1),
setNorms(P+1),
......@@ -180,7 +153,6 @@ namespace xerus { namespace impl_TrASD {
bestTestResidual(std::numeric_limits<double>::max()),
prevRanks(tracking, outX.ranks()),
residuals(tracking, std::numeric_limits<double>::max()),
leftStack(degree, std::vector<Tensor>(numMeasurments)),
rightStack(degree, std::vector<Tensor>(numMeasurments))
......@@ -390,7 +362,7 @@ namespace xerus { namespace impl_TrASD {
}
void finish(const size_t _iteration) {
void finish() {
for(size_t i = 0; i < bestX.degree(); i++) {
if(i == bestX.corePosition) {
outX.set_component(i, bestX.get_average_core());
......@@ -399,7 +371,7 @@ namespace xerus { namespace impl_TrASD {
}
}
LOG(ADF, "Residual decrease from " << std::scientific << 0.0 /* TODO */ << " to " << std::scientific << residuals.back() << " in " << _iteration << " iterations.");
LOG(ASD, "Residual decrease from " << std::scientific << 0.0 /* TODO */ << " to " << std::scientific << current_residual() << " in " << current_iteration() << " iterations.");
}
......@@ -422,19 +394,17 @@ namespace xerus { namespace impl_TrASD {
contract(tmp, setCore, rightStack[corePosition+1][i], 1);
contract(tmp2, measurments.positions[i][corePosition], tmp, 1);
setResiduals[p] += misc::sqr(tmp2[0] - measurments.measuredValues[i]);
// const double test = tmp2[0];
contract(tmp, avgCore, rightStack[corePosition+1][i], 1);
contract(tmp2, measurments.positions[i][corePosition], tmp, 1);
optResidual += misc::sqr(tmp2[0] - measurments.measuredValues[i]);
// LOG(bla, (test-tmp2[0])/test);
REQUIRE(tmp2.size == 1, "IE");
}
setResiduals[p] = std::sqrt(setResiduals[p])/setNorms[p];
}
optResidual = std::sqrt(optResidual)/optNorm;
double testResidual = 0.0;
double testResidual = 0.0;
for(const auto i : sets[P]) {
contract(tmp, avgCore, rightStack[corePosition+1][i], 1);
contract(tmp2, measurments.positions[i][corePosition], tmp, 1);
......@@ -450,7 +420,6 @@ namespace xerus { namespace impl_TrASD {
void solve() {
perfData.start();
size_t nonImprovementCounter = 0;
// Build inital right stack
REQUIRE(x.corePosition == 0, "Expecting core position to be 0.");
......@@ -458,46 +427,42 @@ namespace xerus { namespace impl_TrASD {
update_right_stack(corePosition);
}
for(size_t iteration = 0; maxIterations == 0 || iteration < maxIterations; ++iteration) {
while(!reached_stopping_criteria()) {
double optResidual, testResidual;
std::vector<double> setResiduals;
std::tie(optResidual, testResidual, setResiduals) = calc_residuals();
residuals.push_back(optResidual);
// prevRanks.push_back(x.ranks());
if(P == 1 || testResidual < 0.99*bestTestResidual) {
if(P == 1) {
bestX = x;
bestTestResidual = optResidual;
} else if(testResidual < bestTestResidual) {
bestX = x;
bestTestResidual = testResidual;
nonImprovementCounter = 0;
} else {
nonImprovementCounter++;
}
make_step(bestTestResidual);
// prevRanks.push_back(x.ranks());
// LOG(ADFx, "Residual " << std::scientific << residuals.back() << " " << /*setResiduals*/ -1 << ". NonImpCnt: " << nonImprovementCounter << ", Controlset: " << testResidual << ". Ranks: " << x.ranks() << ". DOFs: " << x.dofs() << ". Norm: " << frob_norm(x.get_average_core()));
// LOG(ASD, "Residual " << std::scientific << optResidual << " " << /*setResiduals*/ -1 << ". Controlset: " << testResidual << ". Ranks: " << x.ranks() << ". DOFs: " << x.dofs() << ". Norm: " << frob_norm(x.get_average_core()));
// bool maxRankReached = true;
// for(size_t k = 0; k+1 < x.degree(); ++k ) {
// maxRankReached = maxRankReached && (x.rank(k) == maxRanks[k]);
// }
if(P > 1 && nonImprovementCounter > 2) {
rankEps = std::min(0.32, 2*rankEps);
if(reached_convergence_criteria()) {
if(misc::hard_equal(rankEps, maxRankEps)) { break; }
reset_convergence_buffer();
rankEps = std::min(maxRankEps, 2*rankEps);
perfData << rankEps;
}
if(optResidual < targetRelativeResidual || nonImprovementCounter >= 25 || ( P == 1 && residuals.back() > std::pow(minimalResidualNormDecrease, tracking)*residuals[0])) {
finish(iteration);
return; // We are done!
}
perfData.add(iteration, std::vector<double>{optResidual, testResidual} | setResiduals, x.get_average_tt(), 0);
perfData.add(current_iteration(), std::vector<double>{optResidual, testResidual} | setResiduals, x.get_average_tt(), 0);
if(P>1) { shuffle_sets(); }
// Forward sweep
for(size_t corePosition = 0; corePosition+1 < degree; ++corePosition) {
update_core(corePosition);
......@@ -509,7 +474,6 @@ namespace xerus { namespace impl_TrASD {
update_core(degree-1);
// Backward sweep
for(size_t corePosition = degree-1; corePosition > 0; --corePosition) {
update_core(corePosition);
......@@ -521,7 +485,7 @@ namespace xerus { namespace impl_TrASD {
update_core(0);
}
finish(maxIterations);
finish();
}
};
......@@ -529,15 +493,17 @@ namespace xerus { namespace impl_TrASD {
void ASDVariant::operator()(TTTensor& _x, const RankOneMeasurementSet& _measurments, PerformanceData& _perfData) const {
impl_TrASD::InternalSolver<1> solver(_x, _measurments, *this, initialRankEps, _x.ranks(), _perfData);
XERUS_REQUIRE_TEST;
impl_TrASD::InternalSolver<1> solver(*this, _x, _measurments, initialRankEps, _x.ranks(), _perfData);
solver.solve();
}
void ASDVariant::operator()(TTTensor& _x, const RankOneMeasurementSet& _measurments, const std::vector<size_t>& _maxRanks, PerformanceData& _perfData) const {
impl_TrASD::InternalSolver<2> solver(_x, _measurments, *this, initialRankEps, _maxRanks, _perfData);
XERUS_REQUIRE_TEST;
impl_TrASD::InternalSolver<2> solver(*this, _x, _measurments, initialRankEps, _maxRanks, _perfData);
solver.solve();
}
const ASDVariant TRASD(1000, 1e-10, 0.9995);
const ASDVariant TRASD(1000, 1e-10, 0.9999);
} // namespace xerus
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