Commit 3e4d763a authored by Michael Goette's avatar Michael Goette

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

parents 3fff1401 11bdf056
......@@ -8,8 +8,8 @@ CXX = g++
# CXX = clang++
PYTHON2_CONFIG = python2-config
PYTHON3_CONFIG = python3-config
PYTHON2_TEST = pytest2
PYTHON3_TEST = pytest3
PYTEST2 = pytest2
PYTEST3 = pytest3
#=================================================================================================
# C++ Version
......
This diff is collapsed.
......@@ -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);
......
import unittest
import numpy as np
import xerus as xe
def build_TestChop(seed, num_tests, max_order=20, max_dimension=100, max_rank=30):
random = np.random.RandomState(seed)
orders = random.randint(low=1, high=max_order+1, size=num_tests)
dimensions = [random.randint(low=1, high=max_dimension, size=order) for order in orders]
ranks = [random.randint(low=1, high=max_rank, size=(order-1)) for order in orders]
def test(dimensions, ranks):
def test_chop(self):
A = xe.TTTensor.random(dimensions.tolist(), ranks.tolist())
L,l,e,r,R = xe.indices(5)
mirrored = lambda pos: A.degree()-1-pos
for corePosition in range(A.degree()):
Al,Ar = A.chop(corePosition)
Ac = A.get_component(corePosition)
res = xe.TensorNetwork()
res(L^corePosition,e,R^mirrored(corePosition)) << Al(L&1,l) * Ac(l,e,r) * Ar(r,R&1)
# norm(A - res)**2 == norm(A)**2 - 2*inner(A,res) + norm(res)**2
nA = xe.frob_norm(A)
nres = xe.frob_norm(res)
inner = xe.Tensor()
inner() << A(e&0) * res(e&0)
self.assertLessEqual(nA**2 - 2*inner[0] + nres**2, 5e-5)
return test_chop
name = lambda d,r: "test_chop_{}_{}".format(".".join(d.astype(str)), ".".join(r.astype(str)))
odir = {name(d,r): test(d,r) for d,r in zip(dimensions, ranks)}
return type("TestChop", (unittest.TestCase,), odir)
TestChop = build_TestChop(0, 100, 4)
......@@ -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");
......
This diff is collapsed.
......@@ -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