Commit 13648b59 authored by Sebastian Wolf's avatar Sebastian Wolf

Work on measurments and ASD

parent 3a4b9469
Pipeline #904 passed with stages
in 8 minutes and 51 seconds
......@@ -34,13 +34,15 @@ namespace xerus {
class ASDVariant : public OptimizationAlgorithm {
public:
double minRankEps = 1e-8;
double minRankEps = 1e-4;
double epsDecay = 0.8;
// double maxRankEps = 1e-1;
double epsDecay = 1.1;
double controlSetFraction = 0.1;
double initialRankEps = 0.01;
double initialRankEps = 5e-2;
/// Basic constructor
ASDVariant(const size_t _maxIterations, const double _targetRelativeResidual, const double _minimalResidualNormDecrease)
......
......@@ -46,7 +46,7 @@ namespace xerus {
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 = 100;
size_t tracking = 10;
protected:
......
......@@ -105,6 +105,8 @@ namespace xerus {
private:
void create_random_positions(const size_t _numMeasurements, const std::vector<size_t>& _dimensions);
void create_slice_random_positions(const size_t _numMeasurements, const std::vector<size_t>& _dimensions);
};
......
......@@ -34,7 +34,7 @@ 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(largestEntry, 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(ADF, xerus::misc::internal::LOGGING_ON_ERROR)
// XERUS_SET_LOGGING(ADFx, xerus::misc::internal::LOGGING_ON_ERROR)
XERUS_SET_LOGGING(UQ, xerus::misc::internal::LOGGING_ON_ERROR)
/* */
......@@ -113,6 +113,32 @@ namespace xerus { namespace impl_TrASD {
std::vector<std::vector<Tensor>> rightStack;
public:
void shuffle_sets() {
//Reset all sets
for(size_t setId = 0; setId < P; ++setId) {
sets[setId].clear();
}
//Create the P sets used for optimization
std::uniform_int_distribution<size_t> setDist(0, P-1);
for(size_t i = 0; i < numMeasurments; ++i ) {
if(!misc::contains(sets[P], i)) {
sets[setDist(misc::randomEngine)].push_back( i );
}
}
// Calculate the measurment norms
for(size_t setId = 0; setId < P; ++setId) {
setNorms[setId] = 0.0;
for(const auto i : sets[setId]) {
setNorms[setId] += misc::sqr(measurments.measuredValues[i]);
}
setNorms[setId] = std::sqrt(setNorms[setId]);
}
}
InternalSolver( TTTensor& _x,
const RankOneMeasurementSet& _measurments,
const ASDVariant& _optiSettings,
......@@ -137,7 +163,7 @@ namespace xerus { namespace impl_TrASD {
minRankEps(_optiSettings.minRankEps),
epsDecay(_optiSettings.epsDecay),
controlSetFraction(_optiSettings.controlSetFraction),
controlSetFraction(P == 1 ? 0.0 : _optiSettings.controlSetFraction),
perfData(_perfData),
......@@ -159,29 +185,23 @@ namespace xerus { namespace impl_TrASD {
XERUS_REQUIRE(numMeasurments > 0, "Need at very least one measurment.");
XERUS_REQUIRE(measurments.degree() == degree, "Measurment degree must coincide with x degree.");
//Create the P sets used for optimization + one test set
// Create test set
std::uniform_real_distribution<double> stochDist(0.0, 1.0);
std::uniform_int_distribution<size_t> setDist(0, P-1);
optNorm = 0.0;
setNorms[P] = 0.0;
for(size_t i = 0; i < numMeasurments; ++i ) {
if(stochDist(misc::randomEngine) > controlSetFraction) {
sets[setDist(misc::randomEngine)].push_back( i );
} else {
if(stochDist(misc::randomEngine) < controlSetFraction) {
sets[P].push_back(i);
setNorms[P] += misc::sqr(measurments.measuredValues[i]);
} else {
optNorm += misc::sqr(measurments.measuredValues[i]);
}
}
// Calculate the measurment norms
optNorm = 0.0;
for(size_t setId = 0; setId < P+1; ++setId ) {
setNorms[setId] = 0.0;
for(const auto i : sets[setId]) {
setNorms[setId] += misc::sqr(measurments.measuredValues[i]);
}
if(setId < P) { optNorm += setNorms[setId]; }
setNorms[setId] = std::sqrt(setNorms[setId]);
}
optNorm = std::sqrt(optNorm);
setNorms[P] = std::sqrt(setNorms[P]);
shuffle_sets();
}
private:
......@@ -437,10 +457,12 @@ namespace xerus { namespace impl_TrASD {
double optResidual, testResidual;
std::vector<double> setResiduals;
std::tie(optResidual, testResidual, setResiduals) = calc_residuals();
residuals.push_back(optResidual);
prevRanks.push_back(x.ranks());
if(testResidual < 0.9999*bestTestResidual) {
// prevRanks.push_back(x.ranks());
if(P == 1 || testResidual < 0.99*bestTestResidual) {
bestX = x;
bestTestResidual = testResidual;
nonImprovementCounter = 0;
......@@ -448,43 +470,35 @@ namespace xerus { namespace impl_TrASD {
nonImprovementCounter++;
}
// 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()));
if(residuals.back() < targetRelativeResidual || nonImprovementCounter >= 1000) {
finish(iteration);
LOG(adoff, "Target or MaxIter");
return;
// 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.25, 2*rankEps);
perfData << rankEps;
}
if(residuals.back() > std::pow(minimalResidualNormDecrease, tracking)*residuals[0]) {
bool maxRankReached = false;
bool rankMaxed = false;
for(size_t k = 0; k < x.degree()-1; ++k ) {
maxRankReached = maxRankReached || (x.rank( k ) == maxRanks[k]);
rankMaxed = rankMaxed || (x.rank( k ) == prevRanks[0][k]+1);
}
if(misc::hard_equal(rankEps, minRankEps) || maxRankReached) {
LOG(adoff, "MinRankEps or maxRank.");
finish(iteration);
return; // We are done!
}
if(!rankMaxed) {
LOG(ADFx, "Reduce rankEps to " << std::max(minRankEps, epsDecay*rankEps));
rankEps = std::max(minRankEps, epsDecay*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);
if(P>1) { shuffle_sets(); }
// Forward sweep
for(size_t corePosition = 0; corePosition+1 < degree; ++corePosition) {
update_core(corePosition);
x.move_core_right(rankEps, std::min(maxRanks[corePosition], prevRanks[1][corePosition]+1));
x.move_core_right(rankEps, /*std::min(*/maxRanks[corePosition]/*, prevRanks[0][corePosition]+1)*/);
update_left_stack(corePosition);
}
......@@ -493,13 +507,13 @@ namespace xerus { namespace impl_TrASD {
// Backward sweep
for(size_t corePosition = degree-1; corePosition > 0; --corePosition) {
// update_core(corePosition);
update_core(corePosition);
x.move_core_left(rankEps, std::min(maxRanks[corePosition-1], prevRanks[1][corePosition-1]+1));
x.move_core_left(rankEps, /*std::min(*/maxRanks[corePosition-1]/*, prevRanks[0][corePosition-1]+1)*/);
update_right_stack(corePosition);
}
// update_core(0);
update_core(0);
}
finish(maxIterations);
......@@ -520,5 +534,5 @@ namespace xerus { namespace impl_TrASD {
solver.solve();
}
const ASDVariant TRASD(1000, 1e-10, 0.99999);
const ASDVariant TRASD(1000, 1e-10, 0.9995);
} // namespace xerus
......@@ -242,6 +242,8 @@ namespace xerus {
void SinglePointMeasurementSet::create_random_positions(const size_t _numMeasurements, const std::vector<size_t>& _dimensions) {
// create_slice_random_positions(_numMeasurements, _dimensions);
XERUS_REQUIRE(misc::product(_dimensions) >= _numMeasurements, "It's impossible to perform as many measurements as requested. " << _numMeasurements << " > " << _dimensions);
// Create distributions
......@@ -267,6 +269,45 @@ namespace xerus {
}
void SinglePointMeasurementSet::create_slice_random_positions(const size_t _sliceDensity, const std::vector<size_t>& _dimensions) {
// Create distributions
std::vector<std::uniform_int_distribution<size_t>> indexDist;
for (size_t i = 0; i < _dimensions.size(); ++i) {
indexDist.emplace_back(0, _dimensions[i]-1);
}
std::set<std::vector<size_t>, vec_compare> measuredPositions;
std::vector<size_t> multIdx(_dimensions.size());
for(size_t mu = 0; mu < _dimensions.size(); ++mu) {
XERUS_REQUIRE(misc::product(_dimensions) >= _sliceDensity*_dimensions[mu], "It's impossible to perform as many measurements as requested. " << _sliceDensity << " > " << _dimensions);
for(size_t k = 0; k < _dimensions[mu]; ++k) {
size_t added = 0;
while (added < _sliceDensity) {
for (size_t i = 0; i < _dimensions.size(); ++i) {
if(i == mu) {
multIdx[i] = k;
} else {
multIdx[i] = indexDist[i](misc::randomEngine);
}
}
if(!misc::contains(measuredPositions, multIdx)) {
measuredPositions.insert(multIdx);
added++;
}
}
}
}
for(const auto& pos : measuredPositions) {
positions.push_back(pos);
}
measuredValues.resize(positions.size(), 0.0);
}
......@@ -521,7 +562,7 @@ namespace xerus {
void RankOneMeasurementSet::create_random_positions(const size_t _numMeasurements, const std::vector<size_t>& _dimensions) {
using ::xerus::misc::operator<<;
XERUS_REQUIRE(misc::product(_dimensions) >= _numMeasurements, "It's impossible to perform as many measurements as requested. " << _numMeasurements << " > " << _dimensions);
// XERUS_REQUIRE(misc::product(_dimensions) >= _numMeasurements, "It's impossible to perform as many measurements as requested. " << _numMeasurements << " > " << _dimensions);
std::vector<Tensor> randOnePosition(_dimensions.size());
while (positions.size() < _numMeasurements) {
......
......@@ -72,10 +72,10 @@ namespace xerus {
void PerformanceData::reset() {
if (active) {
data.clear();
additionalInformation.clear();
startTime = ~0ul;
stopTime = ~0ul;
data.clear();
additionalInformation.clear();
}
}
......
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