Skip to content

Commit

Permalink
more testing
Browse files Browse the repository at this point in the history
  • Loading branch information
@fbischoff
fbischoff committed Sep 27, 2023
1 parent f436cb8 commit 3809700
Show file tree
Hide file tree
Showing 3 changed files with 199 additions and 53 deletions.
65 changes: 44 additions & 21 deletions src/madness/chem/lowrankfunction.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -274,6 +274,12 @@ struct particle {
for (int i=0; i<PDIM; ++i) dims[i]=p[i];
}

std::string str() const {
std::stringstream ss;
ss << *this;
return ss.str();
}

/// assuming two particles only
bool is_first() const {return dims[0]==0;}
/// assuming two particles only
Expand Down Expand Up @@ -377,24 +383,32 @@ struct LRFunctorF12 : public LRFunctorBase<T,NDIM> {

// \int f(1,2)^2 d1d2 = \int f(1,2)^2 pre(1)^2 post(2)^2 d1 d2
typename Tensor<T>::scalar_type term1 =madness::inner(post*post,f12sq(pre*pre));
return term1;
return sqrt(term1);

}

T operator()(const Vector<double,NDIM>& r) const {

if (f12->info.type==OT_SLATER) {
double gamma=f12->info.mu;
auto split = [](const Vector<double,NDIM>& r) {
Vector<double,LDIM> first, second;
for (int i=0; i<LDIM; ++i) {
first[i]=r[i];
second[i]=r[i+LDIM];
}
double result=a(first)*b(second)*exp(-gamma*(first-second).normf());
return result;
} else {
return 1.0;
}
return std::make_pair(first,second);
};

double gamma=f12->info.mu;
auto [first,second]=split(r);

double result=1.0;
if (a.is_initialized()) result*=a(first);
if (b.is_initialized()) result*=b(first);
if (f12->info.type==OT_SLATER) result*=exp(-gamma*(first-second).normf());
else if (f12->info.type==OT_GAUSS) result*=exp(-gamma* madness::inner(first-second,first-second));
else return 1.0;
return result;

}
};

Expand All @@ -418,9 +432,7 @@ struct LRFunctorPure : public LRFunctorBase<T,NDIM> {
}

typename Tensor<T>::scalar_type norm2() const {
double n=f.norm2();
return n*n;

return f.norm2();
}
};

Expand Down Expand Up @@ -537,10 +549,8 @@ struct LRFunctorPure : public LRFunctorBase<T,NDIM> {
/// l2 norm
typename TensorTypeData<T>::scalar_type norm2() const {
auto tmp1=matrix_inner(world,h,h);
auto gg=copy(g);
compress(world,gg);
auto tmp2=matrix_inner(world,gg,g);
return tmp1.trace(tmp2);
auto tmp2=matrix_inner(world,g,g);
return sqrt(tmp1.trace(tmp2));
}

std::vector<Function<T,LDIM>> get_functions(const particle<LDIM>& p) const {
Expand Down Expand Up @@ -750,7 +760,21 @@ struct LRFunctorPure : public LRFunctorBase<T,NDIM> {
* d = inner(lrf(1,2), lrf(1,2))
*/

/// lrf(1,3) = inner(lrf(1,2), lrf(2,3))
/// lrf(1,3) = inner(full(1,2), lrf(2,3))

/// @param[in] f1 the first function
/// @param[in] f2 the second function
/// @param[in] p1 the integration variable of the first function
/// @param[in] p2 the integration variable of the second function
template<typename T, std::size_t NDIM, std::size_t PDIM>
LowRankFunction<T,NDIM> inner(const Function<T,NDIM>& f1, const LowRankFunction<T,NDIM>& f2,
const particle<PDIM> p1, const particle<PDIM> p2) {
auto result=inner(f2,f1,p2,p1);
std::swap(result.g,result.h);
return result;
}

/// lrf(1,3) = inner(lrf(1,2), full(2,3))

/// @param[in] f1 the first function
/// @param[in] f2 the second function
Expand Down Expand Up @@ -793,7 +817,6 @@ struct LRFunctorPure : public LRFunctorBase<T,NDIM> {
// inner(lrf(1,2) ,lrf(2,3) ) = \sum_ij g1_i(1) <h1_i(2) g2_j(2)> h2_j(3)
auto matrix=matrix_inner(world,f2.get_functions(p2),f1.get_functions(p1));
auto htilde=transform(world,f2.get_functions(p2.complement()),matrix);
print("p1 complement",p1.complement());
auto gg=copy(world,f1.get_functions(p1.complement()));
return LowRankFunction<T,NDIM>(gg,htilde,f1.rank_revealing_tol,f1.orthomethod);
}
Expand All @@ -805,14 +828,14 @@ struct LRFunctorPure : public LRFunctorBase<T,NDIM> {
/// @param[in] p1 the integration variable of the first function
/// @param[in] p2 the integration variable of the second function, dummy variable for consistent notation
template<typename T, std::size_t NDIM, std::size_t PDIM>
std::vector<Function<T,NDIM>> inner(const LowRankFunction<T,NDIM>& f1, const std::vector<Function<T,PDIM>>& vf,
std::vector<Function<T,NDIM-PDIM>> inner(const LowRankFunction<T,NDIM>& f1, const std::vector<Function<T,PDIM>>& vf,
const particle<PDIM> p1, const particle<PDIM> p2=particle<PDIM>::particle1()) {
World& world=f1.world;
static_assert(2*PDIM==NDIM);
MADNESS_CHECK(p2.is_first());

// inner(lrf(1,2), f_k(2) ) = \sum_i g1_i(1) <h1_i(2) f_k(2)>
auto matrix=matrix_inner(world,vf,f1.get_functions(p1));
auto matrix=matrix_inner(world,f1.get_functions(p1),vf);
return transform(world,f1.get_functions(p1.complement()),matrix);
}

Expand All @@ -825,12 +848,12 @@ struct LRFunctorPure : public LRFunctorBase<T,NDIM> {
template<typename T, std::size_t NDIM, std::size_t PDIM>
Function<T,NDIM> inner(const LowRankFunction<T,NDIM>& f1, const Function<T,PDIM>& f2,
const particle<PDIM> p1, const particle<PDIM> p2=particle<PDIM>::particle1()) {
return inner(f1,std::vector<Function<T,PDIM>>({f2}),p1,p2);
return inner(f1,std::vector<Function<T,PDIM>>({f2}),p1,p2)[0];
}

template<typename T, std::size_t NDIM, std::size_t LDIM=NDIM/2>
class LowRankFunctionFactory {
public:
public:

const particle<LDIM> p1=particle<LDIM>::particle1();
const particle<LDIM> p2=particle<LDIM>::particle2();
Expand Down
175 changes: 144 additions & 31 deletions src/madness/chem/test_low_rank_function.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -387,6 +387,8 @@ int test_arithmetic(World& world, LowRankFunctionParameters parameters) {
test_output t1("LowRankFunction::arithmetic in dimension " + std::to_string(NDIM));
t1.set_cout_to_terminal();
double thresh=FunctionDefaults<LDIM>::get_thresh();
double thresh_ndim=FunctionDefaults<LDIM>::get_thresh();
print("thresh ldim/ndim",thresh,thresh_ndim);
Function<double,LDIM> phi=FunctionFactory<double,LDIM>(world)
.functor([](const Vector<double,LDIM>& r){return exp(-4.0*inner(r,r));});

Expand All @@ -397,14 +399,16 @@ int test_arithmetic(World& world, LowRankFunctionParameters parameters) {
functor2.f12.reset(GaussOperatorPtr<LDIM>(world,2.0));
functor2.a=phi;

auto builder= LowRankFunctionFactory<double,NDIM>(parameters).set_radius(8)
auto p1=particle<LDIM>::particle1();
auto p2=particle<LDIM>::particle2();

auto builder= LowRankFunctionFactory<double,NDIM>(parameters).set_radius(4)
.set_volume_element(0.1).set_rank_revealing_tol(1.e-10).set_orthomethod("canonical");
auto lrf1=builder.project(functor1);
auto lrf2=builder.project(functor2);

Vector<double,NDIM> r;
r.fill(0.2);
print("lrf1(r)",lrf1(r));

// addition/subtraction
{
Expand Down Expand Up @@ -438,48 +442,155 @@ int test_arithmetic(World& world, LowRankFunctionParameters parameters) {
t1.checkpoint(fabs(n2-refn2)<thresh,"norm2 computation");
}

// inner

// scalar multiplication
{
std::vector<Function<double,LDIM>> arg(3);
for (int i=0; i<3; ++i) arg[i]=FunctionFactory<double,LDIM>(world)
.functor([&i](const Vector<double,LDIM>& r)
{return exp(-r.normf());});
auto p1=particle<LDIM>::particle1();
auto p2=particle<LDIM>::particle2();
auto l1=2.0*lrf1;
t1.checkpoint(fabs(l1(r)-2.0*lrf1(r))<thresh,"oop-place multiplication");
t1.checkpoint(&l1.get_g().front()!=&lrf1.get_g().front(),"subtraction - deep copy");
l1*=0.5;
t1.checkpoint(fabs(l1(r)-lrf1(r))<thresh,"in-place multiplication");

}

return t1.end();
}

template<std::size_t LDIM>
int test_inner(World& world, LowRankFunctionParameters parameters) {

static_assert(LDIM==1 or LDIM==2);
constexpr std::size_t NDIM = 2 * LDIM;
test_output t1("LowRankFunction::test_inner in dimension " + std::to_string(NDIM));
t1.set_cout_to_terminal();
double thresh=FunctionDefaults<LDIM>::get_thresh();
double thresh_ndim=FunctionDefaults<LDIM>::get_thresh();
print("thresh ldim/ndim",thresh,thresh_ndim);
Function<double,LDIM> phi=FunctionFactory<double,LDIM>(world)
.functor([](const Vector<double,LDIM>& r){return exp(-4.0*inner(r,r));});

LRFunctorF12<double,NDIM> functor1;
// functor1.f12.reset(GaussOperatorPtr<LDIM>(world,1.0));
functor1.f12.reset(SlaterOperatorPtr_ND<LDIM>(world,1.0,1.e-4,thresh));
functor1.a=phi;
LRFunctorF12<double,NDIM> functor2;
// functor2.f12.reset(GaussOperatorPtr<LDIM>(world,2.0));
functor2.f12.reset(SlaterOperatorPtr_ND<LDIM>(world,2.0,1.e-4,thresh));
functor2.a=phi;

auto fullrank1=lrf1.reconstruct();
auto fullrank2=lrf2.reconstruct();
t1.checkpoint(true,"prep inner");
auto p1=particle<LDIM>::particle1();
auto p2=particle<LDIM>::particle2();

auto builder= LowRankFunctionFactory<double,NDIM>(parameters).set_radius(4)
.set_volume_element(0.1).set_rank_revealing_tol(1.e-10).set_orthomethod("canonical");
auto lrf1=builder.project(functor1);
auto lrf2=builder.project(functor2);

// reference numbers: (by mathematica)
// f1(x,y) = exp(-a*x^2) * exp(-(x-y)^2)
// f2(x,y) = exp(-a*x^2) * exp(-g (x-y)^2)
// with a=4, g=2
// int f1(x,y),f2(x,z) dx = inner(f1,f2,0,0) : norm^2 = Pi^2/(2 Sqrt[2] Sqrt[a gamma] Sqrt[1 + 2 a + gamma]) = 0.37197471167788324677
// int f1(x,y),f2(z,x) dx = inner(f1,f2,0,1) : norm^2 = 0.32972034117743393239
// int f1(y,x),f2(x,z) dx = inner(f1,f2,1,0) : norm^2 = 0.26921553123369812300
// int f1(y,x),f2(z,x) dx = inner(f1,f2,1,1) : norm^2 = 0.35613867236025352322

// inner f(1,2) f(2,3)
auto fullrank1=lrf1.reconstruct();
auto fullrank2=lrf2.reconstruct();
t1.checkpoint(true,"prep inner");
{
std::vector<double> reference={ 0.37197471167788324677, 0.32972034117743393239, 0.26921553123369812300, 0.35613867236025352322};
int counter=0;
for (auto p11 : {p1,p2}) {
for (auto p22 : {p1,p2}) {
auto lhs0=inner(fullrank1,fullrank2,p11.get_tuple(),p22.get_tuple());
auto lhs1=inner(lrf1,fullrank2,p11,p22);
auto lhs2=inner(lrf1,lrf2,p11,p22);
double l0=lhs0.norm2();
double l1=lhs1.norm2();
double l2=lhs2.norm2();
double ref=reference[counter];
if (LDIM==2) ref*=ref;

print("inner(lrf,full,",p11,p22,"): ",fullrank1.norm2(),lhs1.norm2(),lhs2.norm2());
// full/full
auto lhs1=inner(fullrank1,fullrank2,p11.get_tuple(),p22.get_tuple());
double l1=lhs1.norm2();
t1.checkpoint(fabs(l1*l1-ref)<thresh,"inner(full,full,"+p11.str()+","+p22.str()+")");
double asymmetric_ref=inner(fullrank1,lhs1);
double asymmetric1=inner(fullrank1,lhs1);
t1.checkpoint(fabs(asymmetric1-asymmetric_ref)<thresh,"asymmetric(full,full,"+p11.str()+","+p22.str()+")");

// low rank/full
auto lhs2=inner(lrf1,fullrank2,p11,p22);
double l2=lhs2.norm2();
t1.checkpoint(fabs(l2*l2-ref)<thresh,"inner(lrf,full,"+p11.str()+","+p22.str()+")");
double asymmetric2=inner(lrf1,lhs2);
t1.checkpoint(fabs(asymmetric2-asymmetric_ref)<thresh,"asymmetric(lrf,full,"+p11.str()+","+p22.str()+")");

// full/low rank
auto lhs3=inner(fullrank1,lrf2,p11,p22);
double l3=lhs3.norm2();
t1.checkpoint(fabs(l3*l3-ref)<thresh,"inner(full,lrf,"+p11.str()+","+p22.str()+")");
double asymmetric3=inner(lrf1,lhs3);
t1.checkpoint(fabs(asymmetric3-asymmetric_ref)<thresh,"asymmetric(full,lrf,"+p11.str()+","+p22.str()+")");


// low rank/low rank
auto lhs4=inner(lrf1,lrf2,p11,p22);
double l4=lhs4.norm2();
t1.checkpoint(fabs(l4*l4-ref)<thresh,"inner(lrf,lrf,"+p11.str()+","+p22.str()+")");
double asymmetric4=inner(lrf1,lhs4);
t1.checkpoint(fabs(asymmetric4-asymmetric_ref)<thresh,"asymmetric(lrf,lrf,"+p11.str()+","+p22.str()+")");

print("result norm",p11,p22,"), reference ",l1*l1,l2*l2,l3*l3,l4*l4,ref);
// l2 norm cannot distinguish between f(1,2) and f(2,1)
print("result asym",p11,p22,")",asymmetric1,asymmetric2,asymmetric3,asymmetric4);

counter++;
}
}

// auto lhs1=inner(functor1,arg,p2,p1);
// auto lhs2=inner(lrf1,arg,p2,p1);
// double error=norm2(world,lhs1-lhs2);
}

// scalar multiplication
// inner f(1,2) g(2)
// this is surprisingly inaccurate, but algorithm is correct, the error can be systematically decreased
{
auto l1=2.0*lrf1;
t1.checkpoint(fabs(l1(r)-2.0*lrf1(r))<thresh,"oop-place multiplication");
t1.checkpoint(&l1.get_g().front()!=&lrf1.get_g().front(),"subtraction - deep copy");
l1*=0.5;
t1.checkpoint(fabs(l1(r)-lrf1(r))<thresh,"in-place multiplication");
thresh=FunctionDefaults<LDIM>::get_thresh()*50.0;
// fresh start
lrf1=builder.project(functor1);
fullrank1=FunctionFactory<double,NDIM>(world).functor(functor1);

}
std::vector<Function<double,LDIM>> arg(3);
for (int i=0; i<3; ++i) arg[i]=FunctionFactory<double,LDIM>(world)
.functor([&i](const Vector<double,LDIM>& r)
{return exp(-r.normf());});

std::vector<Function<double,LDIM>> lhs_full1, lhs_full2,lhs_func1,lhs_func2;
for (auto& a : arg) {
lhs_full1.push_back(inner(fullrank1,a,p1.get_tuple(),p1.get_tuple()));
lhs_full2.push_back(inner(fullrank1,a,p2.get_tuple(),p1.get_tuple()));

lhs_func1.push_back(inner(functor1,a,p1,p1));
lhs_func2.push_back(inner(functor1,a,p2,p1));
}
auto lhs_lrf1=inner(lrf1,arg,p1,p1);
auto lhs_lrf2=inner(lrf1,arg,p2,p1);

double norm_func1=norm2(world,lhs_func1);
double norm_func2=norm2(world,lhs_func2);
double norm_full1=norm2(world,lhs_full1);
double norm_full2=norm2(world,lhs_full2);
double norm_lrf1=norm2(world,lhs_lrf1);
double norm_lrf2=norm2(world,lhs_lrf2);
print("norms 1",norm_func1,norm_full1,norm_lrf1);
print("norms 2",norm_func2,norm_full2,norm_lrf2);

double error1=norm2(world,lhs_full1-lhs_func1);
double error2=norm2(world,lhs_full2-lhs_func2);
double error3=norm2(world,lhs_lrf1 -lhs_func1);
double error4=norm2(world,lhs_lrf2 -lhs_func2);

print("error1/2",error1,error2,error3,error4);
// t1.checkpoint(error1<thresh,"inner(full(1,2),g(1)"); // particularly inaccurate, but we're not testing this
// t1.checkpoint(error2<thresh,"inner(full(1,2),g(2)"); // particularly inaccurate, but we're not testing this
t1.checkpoint(error3<thresh,"inner(lrf(1,2),g(1)");
t1.checkpoint(error4<thresh,"inner(lrf(1,2),g(2)");
}
return t1.end();
}

Expand Down Expand Up @@ -527,7 +638,7 @@ int main(int argc, char **argv) {
startup(world, argc, argv);
commandlineparser parser(argc, argv);
int k = parser.key_exists("k") ? std::atoi(parser.value("k").c_str()) : 6;
double thresh = parser.key_exists("thresh") ? std::stod(parser.value("thresh")) : 1.e-5;
double thresh = parser.key_exists("thresh") ? std::stod(parser.value("thresh")) : 1.e-4;
FunctionDefaults<6>::set_tensor_type(TT_2D);

FunctionDefaults<1>::set_thresh(thresh);
Expand Down Expand Up @@ -570,7 +681,9 @@ int main(int argc, char **argv) {
// isuccess+=test_construction_optimization<1>(world,parameters);
// isuccess+=test_construction_optimization<2>(world,parameters);
isuccess+=test_arithmetic<1>(world,parameters);
// isuccess+=test_arithmetic<2>(world,parameters);
isuccess+=test_arithmetic<2>(world,parameters);
isuccess+=test_inner<1>(world,parameters);
isuccess+=test_inner<2>(world,parameters);

// isuccess+=test_lowrank_function(world,parameters);
// isuccess+=test_Kcommutator(world,parameters);
Expand Down
12 changes: 11 additions & 1 deletion src/madness/mra/operator.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1929,7 +1929,17 @@ namespace madness {
}


/// Factory function generating separated kernel for convolution with exp(-mu*r) in 3D
/// Factory function generating separated kernel for convolution with exp(-mu*r) in 3D
/// Note that the 1/(2mu) factor of SlaterF12Operator is not included, this is just the exponential function
template<std::size_t NDIM>
static inline SeparatedConvolution<double, NDIM>* SlaterOperatorPtr_ND(World& world,
double mu, double lo, double eps,
const BoundaryConditions<NDIM>& bc = FunctionDefaults<NDIM>::get_bc(),
int k = FunctionDefaults<NDIM>::get_k()) {
return new SeparatedConvolution<double,NDIM>(world,OperatorInfo(mu,lo,eps,OT_SLATER),bc,k);
}

/// Factory function generating separated kernel for convolution with exp(-mu*r) in 3D
/// Note that the 1/(2mu) factor of SlaterF12Operator is not included, this is just the exponential function
static inline SeparatedConvolution<double, 3>* SlaterOperatorPtr(World& world,
double mu, double lo, double eps,
Expand Down

0 comments on commit 3809700

Please sign in to comment.