1 /*
2   This file is part of MADNESS.
3 
4   Copyright (C) 2007,2010 Oak Ridge National Laboratory
5 
6   This program is free software; you can redistribute it and/or modify
7   it under the terms of the GNU General Public License as published by
8   the Free Software Foundation; either version 2 of the License, or
9   (at your option) any later version.
10 
11   This program is distributed in the hope that it will be useful,
12   but WITHOUT ANY WARRANTY; without even the implied warranty of
13   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14   GNU General Public License for more details.
15 
16   You should have received a copy of the GNU General Public License
17   along with this program; if not, write to the Free Software
18   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 
20   For more information please contact:
21 
22   Robert J. Harrison
23   Oak Ridge National Laboratory
24   One Bethel Valley Road
25   P.O. Box 2008, MS-6367
26 
27   email: harrisonrj@ornl.gov
28   tel:   865-241-3937
29   fax:   865-572-0680
30 
31 
32   $Id$
33 */
34 
35 /// \file mra/testgconv.cc
36 /// \brief Test convolution with Gaussian * polyn
37 
38 #include <madness/mra/mra.h>
39 #include <madness/mra/operator.h>
40 #include <madness/constants.h>
41 
42 using namespace madness;
43 
44 static const int k = 10;
45 static const double thresh = 1e-6;
46 static const double L = 17;
47 
48 // exp(-r^2) / sqrt(pi) = normalized gaussian
g(const coord_1d & r)49 double g(const coord_1d& r) {
50     static const double fac = 1.0/sqrt(constants::pi);
51     return exp(-r[0]*r[0]) * fac;
52 }
53 
54 
55 // g' = exp(-r^2) / sqrt(pi) = derivative of a normalized gaussian
gprime(const coord_1d & r)56 double gprime(const coord_1d& r) {
57     static const double fac = 1.0/sqrt(constants::pi);
58     return -2.* r[0] * exp(-r[0]*r[0]) * fac;
59 }
60 
61 // conv g()
convg(const coord_1d & r)62 double convg(const coord_1d& r) {
63     static const double fac = 1.0/sqrt(2.0*constants::pi);
64     return exp(-r[0]*r[0]*0.5) * fac;
65 }
66 
67 // sqrt(8)*x*exp(-x^2)
h(const coord_1d & r)68 double h(const coord_1d& r) {
69     static const double fac = sqrt(8.0);
70     return exp(-r[0]*r[0]) * fac * r[0];
71 }
72 
73 // g() conv h() == h() conv g()
gconvh(const coord_1d & r)74 double gconvh(const coord_1d& r) {
75     return exp(-0.5*r[0]*r[0]) * r[0];
76 }
77 
78 // D(conv) (g())
conv_prime_g(const coord_1d & r)79 double conv_prime_g(const coord_1d& r) {
80     static const double fac = 1.0/sqrt(2.0*constants::pi);
81     return -exp(-0.5*r[0]*r[0]) * r[0] *fac;
82 }
83 
84 
test_gconv(World & world)85 int test_gconv(World& world) {
86     coord_1d origin(0.0), lo(-L), hi(L);
87     double width = 2.0*L;
88     int success=0;
89 
90     if (world.rank() == 0) print("Test gconv operation");
91 
92     const real_function_1d f = real_factory_1d(world).f(g);
93     double error=f.trace()-1.0;
94     print("error in integral(g) ", error);
95     if (error>FunctionDefaults<1>::get_thresh()) success++;
96     print("success 0 ", success);
97 
98     // convolve with a normalized Gaussian kernel
99     std::vector< std::shared_ptr< Convolution1D<double> > > ops(1);
100     ops[0].reset(new GaussianConvolution1D<double>(k, width/sqrt(constants::pi),
101             width*width, 0, false));
102     real_convolution_1d op(world, ops);
103 
104     real_function_1d opf = op(f);
105     error=opf.trace()-1.0;
106     print("error in integral(op(g)) ", error);
107     if (error>FunctionDefaults<1>::get_thresh()) success++;
108     print("success 1 ", success);
109 
110     real_function_1d exact = real_factory_1d(world).f(convg);
111     print("norm2(g conv g - exact)", (opf-exact).norm2());
112     error=(opf-exact).norm2();
113     if (error>FunctionDefaults<1>::get_thresh()) success++;
114     print("success 2 ", success);
115 
116     real_function_1d q = real_factory_1d(world).f(h);
117     error=q.trace();
118     print("error in integral(h) ", error);
119     if (error>FunctionDefaults<1>::get_thresh()) success++;
120     print("success 3 ", success);
121 
122     error=q.norm2() - sqrt(sqrt(2.0*constants::pi));
123     print("error in norm2(h)", error);
124     if (error>FunctionDefaults<1>::get_thresh()) success++;
125     print("success 4 ", success);
126 
127     real_function_1d opq = op(q);
128     exact = real_factory_1d(world).f(gconvh);
129     error=(opq-exact).norm2();
130     print("norm2(g conv h - exact)", error);
131     if (error>FunctionDefaults<1>::get_thresh()) success++;
132     print("success 5 ", success);
133 
134 
135 
136     // test the convolution with a derivative Gaussian:
137     // result(y) = \int g'(x-y) f(x) = \int g(x-y) f'(x)
138     // where we use
139     // f(x)      = exp(-x^2) / sqrt(pi)
140     // f'(x)     = -2 x exp(-x^2) / sqrt(pi)
141     // result(y) = -y exp(-y/2) / sqrt(2 pi)
142 
143     // the derivative Gaussian convolution kernel:
144     // note the scaling of the coeffs because the derivative operator brings
145     // down the scaling factor of the exponent
146     ops[0].reset(new GaussianConvolution1D<double>(k, 1.0/sqrt(constants::pi),
147             width*width, 1, false));
148 
149     real_convolution_1d oph(world, ops);
150 
151     // this is the result hardwired
152     const real_function_1d convpg=real_factory_1d(world).f(conv_prime_g);
153 
154     // apply the derivative Gaussian on f
155     opq = oph(f);
156 
157     // apply the Gaussian on the derivative of f
158     const real_function_1d fp=real_factory_1d(world).f(gprime);
159     real_function_1d opfp=op(fp);
160 
161     // the error
162     const double error1=(opq-convpg).norm2();
163     const double error2=(opfp-convpg).norm2();
164     print("norm2(conv' g - exact)", error1);
165     print("norm2(conv g' - exact)", error2);
166     if (error1>FunctionDefaults<1>::get_thresh()) success++;
167     print("success 6a ", success);
168     if (error2>FunctionDefaults<1>::get_thresh()) success++;
169     print("success 6b ", success);
170 
171     plot_line("opf.dat", 1001, lo, hi, q, opq, convpg);
172 
173     world.gop.fence();
174     return success;
175 }
176 
177 
main(int argc,char ** argv)178 int main(int argc, char**argv) {
179     initialize(argc,argv);
180     World world(SafeMPI::COMM_WORLD);
181 
182     int success=0;
183     try {
184         startup(world,argc,argv);
185 
186         FunctionDefaults<1>::set_cubic_cell(-L,L);
187         FunctionDefaults<1>::set_k(k);
188         FunctionDefaults<1>::set_thresh(thresh);
189         FunctionDefaults<1>::set_truncate_mode(1);
190         FunctionDefaults<1>::set_initial_level(5);
191         if (world.rank()==0) {
192         	print(" threshold  ", thresh);
193         	print(" polynomial ", k,"\n");
194         }
195         success+=test_gconv(world);
196 
197     }
198     catch (const SafeMPI::Exception& e) {
199         print(e);
200         error("caught an MPI exception");
201     }
202     catch (const madness::MadnessException& e) {
203         print(e);
204         error("caught a MADNESS exception");
205     }
206     catch (const madness::TensorException& e) {
207         print(e);
208         error("caught a Tensor exception");
209     }
210     catch (char* s) {
211         print(s);
212         error("caught a c-string exception");
213     }
214     catch (const char* s) {
215         print(s);
216         error("caught a c-string exception");
217     }
218     catch (const std::string& s) {
219         print(s);
220         error("caught a string (class) exception");
221     }
222     catch (const std::exception& e) {
223         print(e.what());
224         error("caught an STL exception");
225     }
226     catch (...) {
227         error("caught unhandled exception");
228     }
229 
230     world.gop.fence();
231     finalize();
232 
233     return success;
234 }
235 
236