5#ifndef GEMMI_ASUDATA_HPP_
6#define GEMMI_ASUDATA_HPP_
23 void add_point(std::complex<double> x, std::complex<double> y) {
25 double inv_n = 1.0 /
n;
26 double weight = (
n - 1.0) * inv_n;
27 std::complex<double> dx = x -
mean_x;
28 std::complex<double> dy = y -
mean_y;
29 sum_xx += weight * std::norm(dx);
30 sum_yy += weight * std::norm(dy);
31 sum_xy += weight * (dx * std::conj(dy));
35 void add_point(std::complex<float> x, std::complex<float> y) {
36 add_point(std::complex<double>(x), std::complex<double>(y));
44template<
typename Func,
typename T>
46 const std::vector<T>& b,
50 while (r1 != a.end() && r2 != b.end()) {
51 if (r1->hkl == r2->hkl) {
55 }
else if (r1->hkl < r2->hkl) {
66 const std::vector<T>& b) {
77 const std::vector<T>& b) {
90 if (x.value == y.value)
117void move_to_asu(
const GroupOps&,
const Miller& hkl,
int, HklValue<T>& hkl_value) {
122void move_to_asu(
const GroupOps& gops,
const Miller& hkl,
int isym,
123 HklValue<std::complex<R>>& v) {
126 const Op& op = gops.sym_ops[(isym - 1) / 2];
127 double shift = op.phase_shift(hkl);
129 double phase = std::arg(v.value) + shift;
130 v.value = std::polar(std::abs(v.value), (R)phase);
133 v.value.imag(-v.value.imag());
139 std::vector<HklValue<T>>
v;
144 size_t size()
const {
return v.size(); }
146 double get_f(
size_t n)
const {
return std::abs(
v[n].value); }
147 double get_phi(
size_t n)
const {
return std::arg(
v[n].value); }
151 if (!std::is_sorted(
v.begin(),
v.end()))
152 std::sort(
v.begin(),
v.end());
157 fail(
"AsuData::ensure_asu(): space group not set");
170 template<
typename DataProxy>
173 std::size_t col =
proxy.column_index(label);
174 unit_cell_ =
proxy.unit_cell();
175 spacegroup_ =
proxy.spacegroup();
177 auto num =
proxy.get_num(
i + col);
178 if (!std::isnan(num))
179 v.push_back({
proxy.get_hkl(
i), (
T)num});
188 template<
int N,
typename DataProxy>
191 std::array<std::size_t, N>
cols;
192 for (
int i = 0;
i <
N; ++
i)
194 unit_cell_ =
proxy.unit_cell();
195 spacegroup_ =
proxy.spacegroup();
196 using Val =
typename T::value_type;
198 std::array<Val, N>
nums;
199 for (
int j = 0;
j <
N; ++
j)
201 if (!std::any_of(
nums.begin(),
nums.end(), [](
Val f) { return std::isnan(f); })) {
204 set_value_from_array(
v.back().value,
nums);
215 static void set_value_from_array(
T& val,
const std::array<T,1>&
nums) { val =
nums[0]; }
216 static void set_value_from_array(
T& val,
const T& nums) { val = nums; }
218 static void set_value_from_array(std::complex<R>& val,
const std::array<R,2>& nums) {
219 R theta = (R)
rad(nums[1]);
220 val = {nums[0] * std::cos(theta), nums[0] * std::sin(theta)};
223 static void set_value_from_array(ValueSigma<R>& val,
const std::array<R,2>& nums) {
229template<
typename T,
int N,
typename Data>
236template<
typename T,
typename Data>
MtzDataProxy data_proxy(const Mtz &mtz)
void for_matching_reflections(const std::vector< T > &a, const std::vector< T > &b, const Func &func)
std::array< int, 3 > Miller
A synonym for convenient passing of hkl.
int count_equal_values(const std::vector< T > &a, const std::vector< T > &b)
constexpr double rad(double angle)
void fail(const std::string &msg)
ComplexCorrelation calculate_hkl_complex_correlation(const std::vector< T > &a, const std::vector< T > &b)
Correlation calculate_hkl_value_correlation(const std::vector< T > &a, const std::vector< T > &b)
AsuData< T > make_asu_data(const Data &data, const std::array< std::string, N > &labels, bool as_is=false)
double get_phi(size_t n) const
const SpaceGroup * spacegroup() const
Miller get_hkl(size_t n) const
void load_values(const DataProxy &proxy, const std::string &label, bool as_is=false)
double get_f(size_t n) const
std::vector< HklValue< T > > v
void ensure_asu(bool tnt_asu=false)
const UnitCell & unit_cell() const
void load_values(const DataProxy &proxy, const std::array< std::string, N > &labels, bool as_is=false)
std::complex< double > mean_y
void add_point(std::complex< double > x, std::complex< double > y)
double mean_ratio() const
std::complex< double > coefficient() const
std::complex< double > mean_x
void add_point(std::complex< float > x, std::complex< float > y)
std::complex< double > sum_xy
void add_point(double x, double y)
bool operator<(const Miller &m) const
bool operator<(const HklValue &o) const
bool is_in(const Op::Miller &hkl) const
std::pair< Op::Miller, int > to_asu(const Op::Miller &hkl, const GroupOps &gops) const
Returns hkl in asu and MTZ ISYM - 2*n-1 for reflections in the positive asu (I+ of a Friedel pair),...
GroupOps operations() const
bool operator==(const ValueSigma &o) const