Commit 946169c4 authored by Philipp  Trunschke's avatar Philipp Trunschke

Merge branch 'development' of git.hemio.de:xerus/xerus into development

parents 3e35960d ec260050
Pipeline #1057 passed with stages
in 7 minutes and 48 seconds
......@@ -34,6 +34,8 @@ namespace xerus { namespace misc {
std::set<std::string> get_files(const std::string& _path);
bool file_exists(const std::string& _path);
bool file_is_empty(const std::string& _filename);
std::string read_file(const std::string& _path);
......
......@@ -218,6 +218,68 @@ static misc::UnitTest tensor_constructors("Tensor", "Constructors", [](){
static misc::UnitTest tensor_solve_failtests("Tensor", "solve_failtests", []() {
xerus::Tensor A = xerus::Tensor::identity({5, 5});
for(size_t k = 0; k < 2; ++k) { // Sparse and dense A
xerus::Tensor x = xerus::Tensor({4});
xerus::Tensor b = xerus::Tensor::random({3});
FAILTEST(solve(x, A, b););
x = xerus::Tensor({5});
b = xerus::Tensor::random({3});
FAILTEST(solve(x, A, b););
x = xerus::Tensor({6});
b = xerus::Tensor::random({6});
FAILTEST(solve(x, A, b););
x = xerus::Tensor({3});
b = xerus::Tensor::random({5});
solve(x, A, b);
TEST(approx_equal(x, b, 1e-10));
x = xerus::Tensor({5});
b = xerus::Tensor::random({5});
solve(x, A, b);
TEST(approx_equal(x, b, 1e-10));
x = xerus::Tensor({4,3});
b = xerus::Tensor::random({3});
FAILTEST(solve(x, A, b););
if(k==1) { // Currently extra orders only works in dense.
x = xerus::Tensor({4,3});
b = xerus::Tensor::random({3,4});
FAILTEST(solve(x, A, b););
x = xerus::Tensor({5,3});
b = xerus::Tensor::random({3,3});
FAILTEST(solve(x, A, b, 1););
x = xerus::Tensor({5,3});
b = xerus::Tensor::random({3,5});
FAILTEST(solve(x, A, b, 1););
x = xerus::Tensor::random({3,7,2});
b = xerus::Tensor::random({5,3});
solve(x, A, b, 1);
TEST(approx_equal(x, b, 1e-10));
x = xerus::Tensor::random({5,5,2});
b = xerus::Tensor::random({5,3,4});
solve(x, A, b, 2);
TEST(approx_equal(x, b, 1e-10));
}
A.use_dense_representation();
}
return 0;
});
static misc::UnitTest tensor_sparse_dense("Tensor", "Sparse_Dense_Conversions", [](){
Tensor n({3,3,3,3});
const size_t dim = 100;
......
......@@ -48,7 +48,7 @@ namespace xerus { namespace uq {
Tensor hermite_evaluation(const double _v, const size_t _basisSize) {
Tensor p({_basisSize});
for (unsigned i = 0; i < _basisSize; ++i) {
p[i] = std::sqrt(1/(/*std::sqrt(2*M_PI)**/boost::math::factorial<double>(i)))*boost::math::hermite(i, _v/std::sqrt(2))/std::pow(2.0, i/2.0);
p[i] = boost::math::hermite(i, _v/std::sqrt(2))/std::pow(2.0, i/2.0);
}
return p;
}
......@@ -63,6 +63,22 @@ namespace xerus { namespace uq {
}
Tensor normalization_matrix(const PolynomBasis _basisType, const size_t _basisSize) {
Tensor normMat({_basisSize, _basisSize});
if(_basisType == PolynomBasis::Hermite) {
for(uint i = 0; i < _basisSize; ++i) {
normMat[{i, i}] = boost::math::factorial<double>(i);
}
} else {
for(size_t i = 0; i < _basisSize; ++i) {
normMat[{i, i}] = 1.0/(2*double(i)+1);
}
}
return normMat;
}
std::tuple<Tensor, Tensor, Tensor> det_moments(const TTTensor& _x, const PolynomBasis _basisType) {
REQUIRE(_x.corePosition == 0, "Invalid core position.");
......@@ -72,18 +88,16 @@ namespace xerus { namespace uq {
while(m1TT.degree() > 1) {
m1TT.fix_mode(1, 0);
}
double m1Factor = 1.0;
if(_basisType == PolynomBasis::Hermite) {
// m1Factor = 1/(std::sqrt(2*M_PI));
// m1Factor = 1.0/10.0;
}
Tensor m1 = m1Factor*Tensor(m1TT);
Tensor m1 = Tensor(m1TT);
// M2
Tensor spacialCmp = _x.get_component(0);
spacialCmp.reinterpret_dimensions({spacialCmp.dimensions[1], spacialCmp.dimensions[2]});
const auto sqrExpectation = contract(entrywise_product(spacialCmp, spacialCmp), xerus::Tensor::ones({spacialCmp.dimensions[1]}), 1);
Tensor m2 = sqrExpectation;
const Index l1, l2, c1, c2, r1, r2;
Tensor m2 = Tensor::identity({1,1});
for(size_t k = _x.degree()-1; k > 0; --k) {
m2(l1, l2) = _x.get_component(k)(l1, c1, r1)*normalization_matrix(_basisType, _x.get_component(k).dimensions[1])(c1, c2)*m2(r1, r2)*_x.get_component(k)(l2, c2, r2);
}
m2(c1, c2) = _x.get_component(0)(l1, c1, r1)*m2(r1, r2)*_x.get_component(0)(l1, c2, r2);
m2(l1) = Tensor::kronecker({_x.dimensions[0],_x.dimensions[0],_x.dimensions[0]})(l1, c1, c2)*m2(c1, c2);
// M3
Tensor m3(m2.dimensions);
......
......@@ -45,7 +45,7 @@ namespace xerus { namespace uq { namespace impl_uqRaAdf {
const double targetResidual;
const size_t maxRank = 50;
const double minRankEps = 1e-8;
const double minRankEps = 1e-10;
const double epsDecay = 0.8;
const double convergenceFactor = 0.995;
......
......@@ -50,6 +50,10 @@ namespace xerus { namespace misc {
return files;
}
bool file_exists(const std::string& _path) {
return boost::filesystem::exists(_path);
}
bool file_is_empty(const std::string& _filename) {
std::ifstream pFile(_filename);
REQUIRE(pFile.is_open(), "IE: Failed to open file");
......
......@@ -1591,9 +1591,6 @@ namespace xerus {
const size_t degM = _B.degree() - _extraDegree;
const size_t degN = _A.degree() - degM;
REQUIRE(_A.degree() == degM+degN, "Inconsistent dimensions.");
REQUIRE(_B.degree() == degM+_extraDegree, "Inconsistent dimensions.");
// Make sure X has right dimensions
if( _X.degree() != degN + _extraDegree
|| !std::equal(_X.dimensions.begin(), _X.dimensions.begin() + degN, _A.dimensions.begin() + degM)
......@@ -1605,6 +1602,9 @@ namespace xerus {
_X.reset(std::move(newDimX), Tensor::Representation::Dense, Tensor::Initialisation::None);
}
REQUIRE(std::equal(_A.dimensions.begin(), _A.dimensions.begin() + degM, _B.dimensions.begin()), "A and b have incompatible dimensions. A: " << _A.dimensions << ", b: " << _B.dimensions);
XERUS_INTERNAL_CHECK(std::equal(_A.dimensions.begin() + degM, _A.dimensions.end(), _X.dimensions.begin()), "A and b have incompatible dimensions. A: " << _A.dimensions << ", b: " << _B.dimensions);
XERUS_INTERNAL_CHECK(std::equal(_B.dimensions.begin() + degM, _B.dimensions.end(), _X.dimensions.begin()+degN), "A and b have incompatible dimensions. A: " << _A.dimensions << ", b: " << _B.dimensions);
// Calculate multDimensions
const size_t m = misc::product(_A.dimensions, 0, degM);
......@@ -1662,14 +1662,13 @@ namespace xerus {
solve_least_squares(_X, _A, _B, _extraDegree);
return;
}
REQUIRE(&_X != &_B && &_X != &_A, "Not supportet yet");
REQUIRE(&_X != &_B && &_X != &_A, "x=b and x=a is not supported yet.");
const size_t degM = _B.degree() - _extraDegree;
const size_t degN = _A.degree() - degM;
REQUIRE(_A.degree() == degM+degN, "Inconsistent dimensions.");
REQUIRE(_B.degree() == degM+_extraDegree, "Inconsistent dimensions.");
// Make sure X has right dimensions
// Make sure X has the right dimensions
if( _X.degree() != degN + _extraDegree
|| !std::equal(_X.dimensions.begin(), _X.dimensions.begin() + degN, _A.dimensions.begin() + degM)
|| !std::equal(_X.dimensions.begin()+ degN, _X.dimensions.end(), _B.dimensions.begin() + degM))
......@@ -1680,13 +1679,20 @@ namespace xerus {
_X.reset(std::move(newDimX), Tensor::Representation::Dense, Tensor::Initialisation::None);
}
REQUIRE(std::equal(_A.dimensions.begin(), _A.dimensions.begin() + degM, _B.dimensions.begin()), "A and b have incompatible dimensions. A: " << _A.dimensions << ", b: " << _B.dimensions);
XERUS_INTERNAL_CHECK(std::equal(_A.dimensions.begin() + degM, _A.dimensions.end(), _X.dimensions.begin()), "A and b have incompatible dimensions. A: " << _A.dimensions << ", b: " << _B.dimensions);
XERUS_INTERNAL_CHECK(std::equal(_B.dimensions.begin() + degM, _B.dimensions.end(), _X.dimensions.begin()+degN), "A and b have incompatible dimensions. A: " << _A.dimensions << ", b: " << _B.dimensions);
// Calculate multDimensions
const size_t m = misc::product(_A.dimensions, 0, degM);
const size_t n = misc::product(_A.dimensions, degM, degM+degN);
const size_t p = misc::product(_B.dimensions, degM, degM+_extraDegree);
REQUIRE(_B.size == m*p, "A and b have incompatible dimensions. A: " << _A.dimensions << ", b: " << _B.dimensions);
XERUS_INTERNAL_CHECK(_X.size == n*p, "Invalid dimension of x");
// Note that A isdense here
// Note that A is dense here
if(_B.is_dense()) {
blasWrapper::solve(
_X.override_dense_data(),
......
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