Skip to content

Commit

Permalink
separating lowrankfunction and ccpairfunction
Browse files Browse the repository at this point in the history
  • Loading branch information
@fbischoff
fbischoff committed Sep 5, 2023
1 parent e545a63 commit 2dab755
Show file tree
Hide file tree
Showing 3 changed files with 134 additions and 1,080 deletions.
1 change: 1 addition & 0 deletions src/madness/chem/lowrankfunction.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,7 @@


#include<madness/mra/mra.h>
#include<madness/world/timing_utilities.h>
#include<madness/chem/electronic_correlation_factor.h>
#include <random>

Expand Down
316 changes: 19 additions & 297 deletions src/madness/chem/test_ccpairfunction.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,35 +16,6 @@

using namespace madness;

struct XParameters : QCCalculationParametersBase {

XParameters() : QCCalculationParametersBase() {

// initialize with: key, value, comment (optional), allowed values (optional)
initialize<double>("radius",2.0,"the radius");
initialize<double>("volume_element",0.1,"volume covered by each grid point");
initialize<long>("rank",500,"the number of grid points in random grids");
initialize<bool>("stable_power_iteration",true,"use stable power iteration algorithm (orthonormalize)");
initialize<double>("tol",1.e-5,"rank-reduced choleski tolerance");
initialize<std::string>("gridtype","random","the grid type",{"random","cartesian"});
initialize<std::string>("rhsfunctiontype","exponential","the type of function",{"delta","exponential"});
initialize<int>("optimize",0,"number of optimization iterations");
}

void read_and_set_derived_values(World& world, const commandlineparser& parser, std::string tag) {
read_input_and_commandline_options(world,parser,tag);
}

double radius() const {return get<double>("radius");}
double volume_element() const {return get<double>("volume_element");}
double tol() const {return get<double>("tol");}
long rank() const {return get<long>("rank");}
bool stable_power_iteration() const {return get<bool>("stable_power_iteration");}
int optimize() const {return get<int>("optimize");}
std::string gridtype() const {return get<std::string>("gridtype");}
std::string rhsfunctiontype() const {return get<std::string>("rhsfunctiontype");}
};


bool longtest=false;
struct data {
Expand Down Expand Up @@ -161,231 +132,6 @@ class GaussianPotential : public FunctionFunctorInterface<double,3> {
}
};

int test_lowrank_function3(World& world, XParameters& parameters) {
test_output t1("CCPairFunction::low rank function");
t1.set_cout_to_terminal();
madness::default_random_generator.setstate(int(cpu_time())%4149);
madness::default_random_generator.setstate(int(cpu_time())%4149);

constexpr std::size_t LDIM=3;
constexpr std::size_t NDIM=2*LDIM;
print("eps, k, NDIM",FunctionDefaults<NDIM>::get_thresh(),FunctionDefaults<NDIM>::get_k(),NDIM);

parameters.print("grid","end");

Vector<double,LDIM> offset;
offset.fill(0.0);
// Function<double,LDIM> phi1=FunctionFactory<double,LDIM>(world).functor([](const Vector<double,LDIM>& r)
// { return exp(-r.normf());});
Function<double,LDIM> phi1=FunctionFactory<double,LDIM>(world).functor([](const Vector<double,LDIM>& r)
{ return 1.0;});
Function<double,LDIM> phi2=FunctionFactory<double,LDIM>(world).functor([&offset](const Vector<double,LDIM>& r)
{ return exp(-1.0*(r-offset).normf());});
Function<double,LDIM> one=FunctionFactory<double,LDIM>(world)
.functor([](const Vector<double,LDIM>& r) { return 1.0;});

std::shared_ptr<real_convolution_3d> f12(SlaterOperatorPtr(world,1.0,1.e-6,FunctionDefaults<LDIM>::get_thresh()));

// auto f = [](const coord_6d& r) {
// coord_3d r1={r[0],r[1],r[2]};
// coord_3d r2={r[3],r[4],r[5]};
// return exp(-(r1-r2).normf() -r2.normf());
// };

auto compute_result = [&world, &one](const auto& lrf) {
real_function_3d result=real_factory_3d(world);
for (int r=0; r<lrf.rank(); ++r) result+=lrf.g[r]*inner(one,lrf.h[r]);
return result;
};
auto compute_relative_error = [&world,&parameters](const auto reference, const auto result, const auto lrf) {
auto diff=reference-result;
double refnorm=reference.norm2();
double resultnorm=result.norm2();
double error=diff.norm2();
return error/refnorm;
};

auto output=[&parameters] (const double projection_error, const long projection_rank, const double projection_time,
const double optimized_error, const long optimized_rank, const double optimized_time) {
print("error",parameters.radius(),parameters.gridtype(),parameters.rhsfunctiontype(),parameters.volume_element(),
parameters.tol(),
projection_error,projection_rank,projection_time,
optimized_error,optimized_rank,optimized_time);
};

// \phi(1) \bar \phi(1) = \intn phi(1) \phi(2) f(1,2) d2
auto reference = phi1* (*f12)(phi2);
output(0.0,0.0,0.0,0.0,0.0,0.0);

LowRank<double,6> lrf(f12,copy(phi1),copy(phi2));
lrf.stable_power_iteration=parameters.stable_power_iteration();
// plot_plane<6>(world,lrf.lrfunctor,"plot_f12_r2",PlotParameters(world).set_plane({"x1","x4"}));
// lrf.project(parameters.rank(),parameters.radius(),parameters.gridtype(),parameters.rhsfunctiontype());
double cpu0=cpu_time();
lrf.project(parameters.volume_element(),parameters.radius(),parameters.gridtype(),parameters.rhsfunctiontype(),parameters.tol());
double cpu1=cpu_time();

// compare
// \phi(1) \bar \phi(1) = \intn phi(1) \phi(2) f(1,2) d2
// = \int \sum_r g_r(1) h_r(2) d2
// = \sum_r g_r(1) <\phi|h_r>
real_function_3d result=compute_result(lrf);
double projection_error=compute_relative_error(reference,result,lrf);
output(projection_error,lrf.rank(),cpu1-cpu0,0.0,0.0,0.0);

double cpu2=cpu_time();
lrf.optimize(parameters.optimize());
double cpu3=cpu_time();
result=compute_result(lrf);
double optimization_error=compute_relative_error(reference,result,lrf);
output(projection_error,lrf.rank(),cpu1-cpu0,optimization_error,lrf.rank(),cpu3-cpu2);


return t1.end();

}

int test_lowrank_function(World& world) {
test_output t1("CCPairFunction::low rank function");
t1.set_cout_to_terminal();
madness::default_random_generator.setstate(int(cpu_time())%4149);

constexpr std::size_t LDIM=1;
constexpr std::size_t NDIM=2*LDIM;
print("eps, k, NDIM",FunctionDefaults<NDIM>::get_thresh(),FunctionDefaults<NDIM>::get_k(),NDIM);

Function<double,NDIM> f12=FunctionFactory<double,NDIM>(world).functor([&LDIM](const Vector<double,NDIM>& r)
{
Vector<double,LDIM> r1,r2;
for (int i=0; i<LDIM; ++i) {
r1[i]=r[i];
r2[i]=r[i+LDIM];
}
return exp(-(r1-r2).normf());//* exp(-0.2*inner(r1,r1));
});
Function<double,NDIM> phi0=FunctionFactory<double,NDIM>(world).functor([&LDIM](const Vector<double,NDIM>& r)
{
Vector<double,LDIM> r1,r2;
for (int i=0; i<LDIM; ++i) {
r1[i]=r[i];
r2[i]=r[i+LDIM];
}
return exp(-(r1).normf());//* exp(-0.2*inner(r1,r1));
});
Function<double,NDIM> phi1=FunctionFactory<double,NDIM>(world).functor([&LDIM](const Vector<double,NDIM>& r)
{
Vector<double,LDIM> r1,r2;
for (int i=0; i<LDIM; ++i) {
r1[i]=r[i];
r2[i]=r[i+LDIM];
}
return exp(-(r2).normf());//* exp(-0.2*inner(r1,r1));
});
Function<double,NDIM> one=FunctionFactory<double,NDIM>(world)
.functor([&LDIM](const Vector<double,NDIM>& r){ return 1.0; });

// f(1,2) = f12 * phi(1) * phi(2)
Function<double,NDIM> f = phi0*phi1 *f12;
print("lhs = phi0 * phi1; rhs=phi0*f12");
auto lhs=phi0*phi1;
auto rhs=phi0*f12;
auto rhsnorm=rhs.norm2();


phi0.reconstruct();
f12.reconstruct();
Function<double,NDIM> reference=inner(lhs,rhs,{0},{0});
LowRank<double,NDIM> lrf(rhs);
lrf.project(50l,3.0,"random","delta");
lrf.optimize();

double ferror=lrf.error();
print("fnorm, error, rel. error",rhsnorm,ferror,ferror/rhsnorm);


auto result=inner(lhs,lrf,{0},{0});
double refnorm=reference.norm2();

auto reconstruct=result.reconstruct();
auto diff=reference-reconstruct;

double resultnorm=reconstruct.norm2();
double diffnorm=diff.norm2();
print("refnorm, resultnorm ",refnorm, resultnorm);
print("resultnorm, error, rel. error",resultnorm,diffnorm,diffnorm/resultnorm);
print("final - one: k,thresh ",FunctionDefaults<2>::get_k(),FunctionDefaults<2>::get_thresh(),
ferror/rhsnorm, diffnorm/resultnorm);
print("\n");

plot<NDIM>({reference, reconstruct, diff}, "f_and_approx", std::vector<std::string>({"adsf", "asdf", "diff"}));

// return 0;
if (0) {
Function<double,NDIM> f;
double fnorm = f.norm2();
print("norm(2D-f)", fnorm);
t1.checkpoint(true, "hi-dim projection");

long rank = 50;
double radius = 5.0;

LowRank<double, NDIM> lr(copy(f));
lr.project(rank, radius,"random","delta");
double err = lr.error();
print("error in f_approx", err);
lr.orthonormalize(true);
t1.checkpoint(true, "start stuff");

err = lr.error();
print("error in f_approx", err);
t1.checkpoint(true, "compute error");
if (LDIM < 3) {
auto fapprox = lr.reconstruct();
plot<NDIM>({f, fapprox, f - fapprox}, "f_and_approx", std::vector<std::string>({"adsf", "asdf", "diff"}));
t1.checkpoint(true, "plotting");
}

/*
* optimize
*/

lr.optimize(2);
t1.checkpoint(true, "optimization");


err = lr.error();
print("error in f_approx_opt", err);

if (LDIM < 3) {
auto fapprox = lr.reconstruct();
plot<NDIM>({f, fapprox, f - fapprox}, "f_and_approx_opt",
std::vector<std::string>({"adsf", "asdf", "diff"}));
t1.checkpoint(true, "plotting");
}


/// int f(1,2) f12(2,3) d2
f12.print_tree();
auto reference=inner(f,f12,{1},{0});
auto hbar=copy(world,lr.h);
for (auto& hh : hbar) hh=inner(f12,hh,{0},{0});
auto result=lr;
result.h=hbar;
auto reference_approx=result.reconstruct();
result.lrfunctor.f=reference;
double error=result.error();
print("error ",error);
plot<NDIM>({reference, reference_approx, reference-reference_approx}, "contraction", std::vector<std::string>({"adsf", "asdf", "diff"}));
}

// auto result=




return t1.end();
}

int test_constructor(World& world, std::shared_ptr<NuclearCorrelationFactor> ncf, const Molecule& molecule,
const CCParameters& parameter) {
test_output t1("CCPairFunction::constructor");
Expand Down Expand Up @@ -1241,35 +987,12 @@ int main(int argc, char **argv) {
madness::World& world = madness::initialize(argc, 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-4;
FunctionDefaults<6>::set_tensor_type(TT_2D);

FunctionDefaults<1>::set_thresh(thresh);
FunctionDefaults<2>::set_thresh(thresh);
FunctionDefaults<3>::set_thresh(thresh);
FunctionDefaults<4>::set_thresh(thresh);
FunctionDefaults<5>::set_thresh(thresh);
FunctionDefaults<6>::set_thresh(1.e-4);

FunctionDefaults<1>::set_k(k);
FunctionDefaults<2>::set_k(k);
FunctionDefaults<3>::set_k(k);
FunctionDefaults<4>::set_k(k);
FunctionDefaults<5>::set_k(k);
FunctionDefaults<6>::set_k(k);

FunctionDefaults<1>::set_cubic_cell(-10.,10.);
FunctionDefaults<2>::set_cubic_cell(-10.,10.);
FunctionDefaults<3>::set_cubic_cell(-10.,10.);
FunctionDefaults<4>::set_cubic_cell(-10.,10.);
FunctionDefaults<6>::set_cubic_cell(-10.,10.);

FunctionDefaults<2>::set_tensor_type(TT_FULL);
FunctionDefaults<3>::set_thresh(1.e-5);
FunctionDefaults<3>::set_cubic_cell(-1.0,1.0);
FunctionDefaults<6>::set_cubic_cell(-1.0,1.0);
print("numerical parameters: k, eps(3D), eps(6D)", FunctionDefaults<3>::get_k(), FunctionDefaults<3>::get_thresh(),
FunctionDefaults<6>::get_thresh());
XParameters parameters;
parameters.read_and_set_derived_values(world,parser,"grid");
int isuccess=0;
#ifdef USE_GENTENSOR

Expand All @@ -1280,25 +1003,24 @@ int main(int argc, char **argv) {
mol.print();
CCParameters ccparam(world, parser);

// data1=data(world,ccparam);
data1=data(world,ccparam);

std::shared_ptr<NuclearCorrelationFactor> ncf = create_nuclear_correlation_factor(world,
mol, nullptr, std::make_pair("slater", 2.0));

isuccess+=test_lowrank_function3(world,parameters);
// isuccess+=test_constructor(world, ncf, mol, ccparam);
// isuccess+=test_operator_apply(world, ncf, mol, ccparam);
// isuccess+=test_transformations(world, ncf, mol, ccparam);
// isuccess+=test_inner(world, ncf, mol, ccparam);
// isuccess+=test_multiply(world, ncf, mol, ccparam);
// isuccess+=test_multiply_with_f12(world, ncf, mol, ccparam);
// isuccess+=test_swap_particles(world, ncf, mol, ccparam);
// isuccess+=test_scalar_multiplication(world, ncf, mol, ccparam);
// isuccess+=test_partial_inner_3d(world, ncf, mol, ccparam);
// isuccess+=test_partial_inner_6d(world, ncf, mol, ccparam);
// isuccess+=test_projector(world, ncf, mol, ccparam);
// FunctionDefaults<3>::set_cubic_cell(-10,10);
// isuccess+=test_helium(world,ncf,mol,ccparam);
mol, nullptr, std::make_pair("slater", 2.0));

isuccess+=test_constructor(world, ncf, mol, ccparam);
isuccess+=test_operator_apply(world, ncf, mol, ccparam);
isuccess+=test_transformations(world, ncf, mol, ccparam);
isuccess+=test_inner(world, ncf, mol, ccparam);
isuccess+=test_multiply(world, ncf, mol, ccparam);
isuccess+=test_multiply_with_f12(world, ncf, mol, ccparam);
isuccess+=test_swap_particles(world, ncf, mol, ccparam);
isuccess+=test_scalar_multiplication(world, ncf, mol, ccparam);
isuccess+=test_partial_inner_3d(world, ncf, mol, ccparam);
isuccess+=test_partial_inner_6d(world, ncf, mol, ccparam);
isuccess+=test_projector(world, ncf, mol, ccparam);
FunctionDefaults<3>::set_cubic_cell(-10,10);
isuccess+=test_helium(world,ncf,mol,ccparam);
data1.clear();
} catch (std::exception& e) {
madness::print("an error occured");
Expand Down
Loading

0 comments on commit 2dab755

Please sign in to comment.