Commit 605f2846 authored by Sebastian Wolf's avatar Sebastian Wolf

Tidy up UQ ADF

parent 635321fe
Pipeline #678 failed with stages
in 3 minutes and 40 seconds
// Xerus - A General Purpose Tensor Library
// Copyright (C) 2014-2016 Benjamin Huber and Sebastian Wolf.
// Copyright (C) 2014-2017 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
......@@ -23,8 +23,6 @@
*/
#include <xerus/algorithms/uqAdf.h>
#include <xerus/indexedTensorMoveable.h>
#include <xerus/misc/basicArraySupport.h>
#include <xerus/misc/simpleNumerics.h>
......@@ -46,8 +44,6 @@ namespace xerus {
TTTensor& x;
std::vector<std::vector<Tensor>> leftStack; // From corePosition 1 to d-2
std::vector<std::vector<Tensor>> rightStack; // From corePosition 1 to d-1
std::vector<std::vector<Tensor>> leftIsStack;
......@@ -98,7 +94,6 @@ namespace xerus {
positions(create_positions(_x, _randomVariables)),
solutions(_solutions),
x(_x),
leftStack(d, std::vector<Tensor>(N)),
rightStack(d, std::vector<Tensor>(N)),
leftIsStack(d, std::vector<Tensor>(N)),
leftOughtStack(d, std::vector<Tensor>(N))
......@@ -129,19 +124,12 @@ namespace xerus {
if(_corePosition > 1) {
contract(tmp, measCmp, true, leftIsStack[_corePosition-1][j], false, 1);
contract(leftIsStack[_corePosition][j], tmp, false, measCmp, false, 1);
contract(leftIsStack[_corePosition][j], tmp, measCmp, 1);
} else { // _corePosition == 1
contract(leftIsStack[_corePosition][j], measCmp, true, measCmp, false, 1);
}
contract(leftOughtStack[_corePosition][j], leftOughtStack[_corePosition-1][j], measCmp, 1);
// if(_corePosition > 1) {
// contract(leftStack[_corePosition][j], leftStack[_corePosition-1][j], measCmp, 1);
// } else { // _corePosition == 1
// leftStack[_corePosition][j] = measCmp;
// }
}
}
}
......@@ -166,107 +154,51 @@ namespace xerus {
}
}
std::vector<Tensor> calc_residual(const size_t _corePosition) {
// LOG(bug, "Calculate residual @ " << _corePosition);
std::vector<Tensor> residual;
Tensor tmp;
if(_corePosition > 0) {
const Tensor shuffledX = reshuffle(x.get_component(_corePosition), {1, 0, 2});
for(size_t j = 0; j < N; ++j) {
// Current node
contract(tmp, positions[_corePosition][j], shuffledX, 1);
// Right part
if(_corePosition < d-1) {
contract(tmp, tmp, rightStack[_corePosition+1][j], 1);
} else {
contract(tmp, tmp, one, 1);
}
// Left Part
if(_corePosition > 1) {
contract(tmp, leftStack[_corePosition-1][j], tmp, 1);
}
// Spacial Part
contract(tmp, x.get_component(0), tmp, 1);
contract(tmp, one, tmp, 1); // TODO weg
residual.push_back(tmp - solutions[j]);
}
} else { // _corePosition == 0
for(size_t j = 0; j < N; ++j) {
// Spacial Part
contract(tmp, x.get_component(0), rightStack[1][j], 1);
contract(tmp, one, tmp, 1); // TODO weg
residual.push_back(tmp - solutions[j]);
}
}
return residual;
}
Tensor calculate_delta(const size_t _corePosition) {
//calculate_projected_gradient
Tensor delta(x.get_component(_corePosition).dimensions);
Tensor dyadComp;
Tensor dyadComp, tmp;
Tensor tmp;
if(_corePosition > 0) {
const Tensor shuffledX = reshuffle(x.get_component(_corePosition), {1, 0, 2});
for(size_t j = 0; j < N; ++j) {
// Calculate common "dyadic part"
Tensor dyadicPart;
if(_corePosition < d-1) {
contract(dyadicPart, positions[_corePosition][j], rightStack[_corePosition+1][j], 0);
} else {
dyadicPart = positions[_corePosition][j];
dyadicPart.reinterpret_dimensions({dyadicPart.dimensions[0], 1}); // Add dangling 1-mode
}
// Calculate "is"
contract(tmp, positions[_corePosition][j], shuffledX, 1);
Tensor isPart;
contract(isPart, positions[_corePosition][j], shuffledX, 1);
if(_corePosition < d-1) {
contract(tmp, tmp, rightStack[_corePosition+1][j], 1);
contract(isPart, isPart, rightStack[_corePosition+1][j], 1);
} else {
tmp.reinterpret_dimensions({tmp.dimensions[0]});
isPart.reinterpret_dimensions({isPart.dimensions[0]});
}
if(_corePosition > 1) { // NOTE: For _corePosition == 1 leftIsStack is the identity
contract(tmp, leftIsStack[_corePosition-1][j], tmp, 1);
contract(isPart, leftIsStack[_corePosition-1][j], isPart, 1);
}
contract(tmp, tmp, positions[_corePosition][j], 0);
if(_corePosition < d-1) {
contract(dyadComp, tmp, rightStack[_corePosition+1][j], 0);
} else {
dyadComp = tmp;
dyadComp.reinterpret_dimensions({dyadComp.dimensions[0], dyadComp.dimensions[1], 1});
}
// Combine with ought part
contract(dyadComp, isPart - leftOughtStack[_corePosition-1][j], dyadicPart, 0);
delta += dyadComp;
// Calculate "Ought"
contract(tmp, leftOughtStack[_corePosition-1][j], positions[_corePosition][j], 0);
if(_corePosition < d-1) {
contract(dyadComp, tmp, rightStack[_corePosition+1][j], 0);
} else {
dyadComp = tmp;
dyadComp.reinterpret_dimensions({dyadComp.dimensions[0], dyadComp.dimensions[1], 1});
}
delta -= dyadComp;
}
} else { // _corePosition == 0
Tensor shuffledX = x.get_component(0);
shuffledX.reinterpret_dimensions({shuffledX.dimensions[1], shuffledX.dimensions[2]});
for(size_t j = 0; j < N; ++j) {
contract(dyadComp, x.get_component(0), rightStack[_corePosition+1][j], 1);
dyadComp.reinterpret_dimensions({x.dimensions[0]});
dyadComp -= solutions[j];
dyadComp.reinterpret_dimensions({1, x.dimensions[0]});
contract(dyadComp, dyadComp, rightStack[_corePosition+1][j], 0);
contract(dyadComp, shuffledX, rightStack[_corePosition+1][j], 1);
contract(dyadComp, dyadComp - solutions[j], rightStack[_corePosition+1][j], 0);
dyadComp.reinterpret_dimensions({1, dyadComp.dimensions[0], dyadComp.dimensions[1]});
delta += dyadComp;
}
}
......@@ -316,26 +248,12 @@ namespace xerus {
return std::sqrt(norm);
}
void update_x(const Tensor& _delta, const double _normAProjGrad, const size_t _corePosition) {
const value_t PyR = misc::sqr(frob_norm(_delta));
x.component(_corePosition) -= (PyR/misc::sqr(_normAProjGrad))*_delta;
}
double residual_norm(const std::vector<Tensor> _resiudal) {
double norm = 0;
for(const auto& r : _resiudal) {
norm += misc::sqr(frob_norm(r));
}
return std::sqrt(norm);
}
double calc_residual_norm(const size_t _corePosition) {
REQUIRE(_corePosition == 0, "Invalid corePosition");
......@@ -352,17 +270,16 @@ namespace xerus {
return std::sqrt(norm);
}
void solve() {
const size_t maxIterations = 10;
for(size_t iteration = 0; maxIterations == 0 || iteration < maxIterations; ++iteration) {
x.require_correct_format();
LOG(start, x.dimensions << " and " << x.ranks());
for(size_t corePosition = 0; corePosition < x.degree(); ++corePosition) {
REQUIRE(x.get_component(corePosition).all_entries_valid(), "Invalid entries at " << corePosition);
}
// x.require_correct_format();
// for(size_t corePosition = 0; corePosition < x.degree(); ++corePosition) {
// REQUIRE(x.get_component(corePosition).all_entries_valid(), "Invalid entries at " << corePosition);
// }
x.move_core(0, true);
// Rebuild right stack
......@@ -372,18 +289,12 @@ namespace xerus {
for(size_t corePosition = 0; corePosition < x.degree(); ++corePosition) {
// const auto residual = calc_residual(corePosition);
if(corePosition == 0) {
// LOG(bla, "Residual @ " << corePosition << " norm: " << std::scientific << residual_norm(residual)/solutionsNorm);
LOG(bla, "Residual norm: " << std::scientific << calc_residual_norm(0)/solutionsNorm);
}
const auto delta = calculate_delta(corePosition);
const auto normAProjGrad = calculate_norm_A_projGrad(delta, corePosition);
update_x(delta, normAProjGrad, corePosition);
// If we have not yet reached the end of the sweep we need to take care of the core and update our stacks
......
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