1 ///////////////////////////////////////////////////////////////////////////////////////////////////
2 /// \file cqmc/numeric/numeric.cpp
3 ///
4 /// \brief implementation file for miscellaneous functions related to numbers
5 ///
6 ///////////////////////////////////////////////////////////////////////////////////////////////////
7
8 //#include <mpi.h>
9 #include "formic/utils/mpi_interface.h"
10
11 #include "formic/utils/zero_one.h"
12 #include "formic/utils/numeric.h"
13 #include "engine_numeric.h"
14
15
16 ///////////////////////////////////////////////////////////////////////////////////////////////////
17 /// \brief computes the unbiased estimate of a ratio of means: <f>_p / <g>_p in which the
18 /// numerator and denominator values are sampled from the same probability distribution p
19 ///
20 /// \param[in] n the number of samples
21 /// \param[in] p the probability weight for each sample
22 /// \param[in] f the numerator samples
23 /// \param[in] g the denominator samples
24 /// \param[out] r on exit, the estimate of the ratio <f>_p / <g>_p
25 /// \param[out] v on exit, the estimate of the variance in the ratio
26 ///
27 ///////////////////////////////////////////////////////////////////////////////////////////////////
unbiased_ratio_of_means(const int n,const double * const p,const S * const f,const T * const g,const bool correct,double & r,double & v)28 template<typename S, typename T> void cqmc::unbiased_ratio_of_means(const int n, const double * const p, const S * const f, const T * const g, const bool correct, double & r, double & v) {
29
30 // compute the normalization, the numerator and denominator means, the means of the squares, and the mean of the products
31 double nm = 0.0; // normalization constant
32 double mf = 0.0; // mean of numerator
33 double mg = 0.0; // mean of denominator
34 double sf = 0.0; // mean of the square of the numerator terms
35 double sfnc = 0.0; // mean of the square of the numerator used in exact sampling
36 double sg = 0.0; // mean of the square of the denominator terms
37 double mp = 0.0; // mean of the product of numerator times denominator
38 for (int i = 0; i < n; i++) {
39 nm += formic::real(p[i]);
40 S x = p[i] * f[i];
41 mf += formic::real(x);
42 if ( correct ) {
43 sf += formic::real(x * f[i]);
44 mp += formic::real(x * g[i]);
45 }
46 else
47 sfnc += formic::real(p[i] * (f[i] / g[i]) * (f[i] / g[i]) * g[i]);
48 x = p[i] * g[i];
49 mg += formic::real(x);
50 if ( correct )
51 sg += formic::real(x * g[i]);
52 }
53 mf /= nm;
54 mg /= nm;
55 sf /= nm;
56 sfnc /= nm;
57 sg /= nm;
58 mp /= nm;
59
60 // compute the numerator and denominator variances and the covariance
61 const double vf = ( sf - mf * mf ) * static_cast<double>(n) / static_cast<double>(n-1);
62 const double vg = ( sg - mg * mg ) * static_cast<double>(n) / static_cast<double>(n-1);
63 const double cv = ( mp - mf * mg ) * static_cast<double>(n) / static_cast<double>(n-1);
64
65 double z = ( correct ? 1.0 : 0.0);
66
67 // compute the unbiased estimate of the ratio of means
68 r = ( mf / mg ) / ( 1.0 + z * ( vg / mg / mg - cv / mf / mg ) / static_cast<double>(n) );
69
70 // compute the unbiased estimate of the variance of the ratio of means
71 if ( correct )
72 v = ( mf * mf / mg / mg / static_cast<double>(n) ) * ( vf / mf / mf + vg / mg / mg - 2.0 * cv / mf / mg );
73 else
74 v = sfnc - mf * mf;
75
76 }
77
78 template void cqmc::unbiased_ratio_of_means(const int n, const double * const p, const double * const f, const double * const g, const bool correct, double & r, double & v);
79 template void cqmc::unbiased_ratio_of_means(const int n,
80 const double * const p,
81 const std::complex<double> * const f,
82 const std::complex<double> * const g,
83 const bool correct,
84 double & r,
85 double & v);
86 template void cqmc::unbiased_ratio_of_means(const int n,
87 const double * const p,
88 const std::complex<double> * const f,
89 const double * const g,
90 const bool correct,
91 double & r,
92 double & v);
93
94
95 ///////////////////////////////////////////////////////////////////////////////////////////////////
96 /// \brief computes the unbiased estimate of a ratio of means: <f>_p / <g>_p in which the
97 /// numerator and denominator values are sampled from the same probability distribution p
98 /// and samples are combined across all processors
99 ///
100 /// \param[in] n the number of samples on this process
101 /// \param[in] p the probability weight for each sample
102 /// \param[in] f the numerator samples
103 /// \param[in] g the denominator samples
104 /// \param[out] r on exit, the estimate of the ratio <f>_p / <g>_p
105 /// \param[out] v on exit, the estimate of the variance in the ratio
106 ///
107 ///////////////////////////////////////////////////////////////////////////////////////////////////
mpi_unbiased_ratio_of_means(const int n,const double * const p,const S * const f,const T * const g,const bool correct,S & r,double & v)108 template<typename S, typename T> void cqmc::mpi_unbiased_ratio_of_means(const int n, const double * const p, const S * const f, const T * const g, const bool correct, S & r, double & v) {
109
110 // compute the normalization, the numerator and denominator means, the means of the squares, and the mean of the products
111 S y[8];
112 y[0] = 0.0; // normalization constant
113 y[1] = 0.0; // mean of numerator
114 y[2] = 0.0; // mean of denominator
115 y[3] = 0.0; // mean of the square of the numerator terms
116 y[4] = 0.0; // mean of the square of the denominator terms
117 y[5] = 0.0; // mean of the product of numerator times denominator
118 y[6] = static_cast<S>(n); // number of samples
119 y[7] = 0.0; // mean of square of numerator terms used in exact sampling(when correct is false)
120
121 // normalization constant
122 double norm = 0.0;
123
124 // accumulate data
125 for (int i = 0; i < n; i++) {
126 norm += p[i];
127 S x = p[i] * f[i];
128 y[1] += x;
129 if ( correct ) {
130 y[3] += x * formic::conj(f[i]);
131 y[5] += x * g[i];
132 }
133 else
134 y[7] += p[i] * f[i] * formic::conj(f[i]);
135 x = p[i] * g[i];
136 y[2] += x;
137 if ( correct )
138 y[4] += x * formic::conj(g[i]);
139 }
140 S z[8];
141 double total_norm = 0.0;
142 formic::mpi::allreduce(&y[0], &z[0], 8, MPI_SUM);
143 formic::mpi::allreduce(&norm, &total_norm, 1, MPI_SUM);
144
145 const S mf = z[1] / total_norm; // mean of numerator
146 const S mg = z[2] / total_norm; // mean of denominator
147 const S sf = z[3] / total_norm; // mean of the square of the numerator terms
148 const S sg = z[4] / total_norm; // mean of the square of the denominator terms
149 const S mp = z[5] / total_norm; // mean of the product of numerator times denominator
150 const S sfnc = z[7] / total_norm; // mean of square of the numerator used in exact sampling
151 const S ns = z[6]; // number of samples
152
153 // compute the numerator and denominator variances and the covariance
154 const S vf = ( sf - mf * mf ) * ns / ( ns - 1.0 );
155 const S vg = ( sg - mg * mg ) * ns / ( ns - 1.0 );
156 const S cv = ( mp - mf * mg ) * ns / ( ns - 1.0 );
157
158 double c = (correct ? 1.0 : 0.0);
159
160 // compute the unbiased estimate of the ratio of means
161 r = ( mf / mg ) / ( 1.0 + c * ( vg / mg / mg - cv / mf / mg ) / ns );
162 //std::cout << vg / mg / mg << " " << cv / mf / mg << std::endl;
163
164 // compute the unbiased estimate of the variance of the ratio of means
165 if ( correct )
166 v = formic::real(( mf * formic::conj(mf) / mg / mg ) * ( vf / mf / formic::conj(mf) + vg / mg / mg - 2.0 * cv / mf / mg ));
167 else
168 v = formic::real(sfnc - mf * formic::conj(mf));
169
170 }
171
172 template void cqmc::mpi_unbiased_ratio_of_means(const int n, const double * const p, const double * const f, const double * const g, const bool correct, double & r, double & v);
173 template void cqmc::mpi_unbiased_ratio_of_means(const int n,
174 const double * const p,
175 const std::complex<double> * const f,
176 const std::complex<double> * const g,
177 const bool correct,
178 std::complex<double> & r,
179 double & v);
180 template void cqmc::mpi_unbiased_ratio_of_means(const int n,
181 const double * const p,
182 const std::complex<double> * const f,
183 const double * const g,
184 const bool correct,
185 std::complex<double> & r,
186 double & v);
187
188
189 ///////////////////////////////////////////////////////////////////////////////////////////////////
190 /// \brief rounds a double to the nearest integer
191 ///
192 /// \param[in] d the number to round
193 ///
194 /// \return the integer nearst to d
195 ///////////////////////////////////////////////////////////////////////////////////////////////////
my_round(const double d)196 int cqmc::my_round(const double d) { return ( d >= 0 ? int(d + 0.5) : int(d - 0.5) ); }
197