/ compute the L2 norm of two complex-valued HDF5 datasets, */ after (a) normalizing each data set by its maximum amplitude*/ and (b) compensating for a constant overall phase factor */ between the datasets. */ /
| 141 | /* between the datasets. */ |
| 142 | /***************************************************************/ |
| 143 | double compare_complex_hdf5_datasets(const char *file1, const char *name1, const char *file2, |
| 144 | const char *name2, int expected_rank = 2, |
| 145 | double *max_dft = 0) { |
| 146 | char dataname[100]; |
| 147 | |
| 148 | // read dataset 1 |
| 149 | h5file f1(file1, h5file::READONLY, false); |
| 150 | int rank1; |
| 151 | size_t *dims1 = new size_t[expected_rank]; |
| 152 | snprintf(dataname, 100, "%s.r", name1); |
| 153 | double *rdata1 = |
| 154 | (double *)f1.read(dataname, &rank1, dims1, expected_rank, false /* single_precision */); |
| 155 | snprintf(dataname, 100, "%s.i", name1); |
| 156 | double *idata1 = |
| 157 | (double *)f1.read(dataname, &rank1, dims1, expected_rank, false /* single_precision */); |
| 158 | if (!rdata1 || !idata1) return -1.0; |
| 159 | |
| 160 | // read dataset 2 |
| 161 | h5file f2(file2, h5file::READONLY, false); |
| 162 | int rank2; |
| 163 | size_t *dims2 = new size_t[expected_rank]; |
| 164 | snprintf(dataname, 100, "%s.r", name2); |
| 165 | double *rdata2 = |
| 166 | (double *)f2.read(dataname, &rank2, dims2, expected_rank, false /* single_precision */); |
| 167 | snprintf(dataname, 100, "%s.i", name2); |
| 168 | double *idata2 = |
| 169 | (double *)f2.read(dataname, &rank2, dims2, expected_rank, false /* single_precision */); |
| 170 | if (!rdata2 || !idata2) return -1.0; |
| 171 | |
| 172 | // check same size |
| 173 | bool same_size = (rank1 == rank2); |
| 174 | for (int d = 0; same_size && d < rank1; d++) |
| 175 | if (dims1[d] != dims2[d]) same_size = false; |
| 176 | if (!same_size) return -1.0; |
| 177 | |
| 178 | // first pass to normalize each dataset to its maximum absolute magnitude; |
| 179 | // we also note the phase difference between the datasets at their points |
| 180 | // of maximum magnitude so we can compensate for this in the comparison below. |
| 181 | size_t length = dims1[0]; |
| 182 | for (int d = 1; d < rank1; d++) |
| 183 | length *= dims1[d]; |
| 184 | |
| 185 | double max_abs1 = 0.0, max_abs2 = 0.0; |
| 186 | double max_arg1 = 0.0, max_arg2 = 0.0; |
| 187 | for (size_t n = 0; n < length; n++) { |
| 188 | std::complex<double> z1 = std::complex<double>(rdata1[n], idata1[n]); |
| 189 | if (abs(z1) > max_abs1) { |
| 190 | max_abs1 = abs(z1); |
| 191 | max_arg1 = arg(z1); |
| 192 | } |
| 193 | std::complex<double> z2 = std::complex<double>(rdata2[n], idata2[n]); |
| 194 | if (abs(z2) > max_abs2) { |
| 195 | max_abs2 = abs(z2); |
| 196 | max_arg2 = arg(z2); |
| 197 | } |
| 198 | } |
| 199 | *max_dft = max_abs1; |
| 200 |
no test coverage detected