00001 #ifndef CCTBX_CRYSTAL_PAIR_TABLES_H
00002 #define CCTBX_CRYSTAL_PAIR_TABLES_H
00003
00004 #include <cctbx/crystal/neighbors_fast.h>
00005 #include <boost/scoped_array.hpp>
00006 #include <map>
00007 #include <set>
00008
00009 namespace cctbx { namespace crystal {
00010
00012 typedef sgtbx::rt_mx pair_sym_op;
00014 typedef std::vector<sgtbx::rt_mx> pair_sym_ops;
00016 typedef std::map<unsigned, pair_sym_ops> pair_sym_dict;
00018 typedef af::shared<pair_sym_dict> pair_sym_table;
00019
00023 inline
00024 af::shared<double>
00025 get_distances(
00026 af::const_ref<pair_sym_dict> const& pair_sym_table,
00027 scitbx::mat3<double> const& orthogonalization_matrix,
00028 af::const_ref<scitbx::vec3<double> > const& sites_frac)
00029 {
00030 CCTBX_ASSERT(sites_frac.size() == pair_sym_table.size());
00031 af::shared<double> distances;
00032 for(unsigned i_seq=0;i_seq<pair_sym_table.size();i_seq++) {
00033 crystal::pair_sym_dict const& pair_sym_dict = pair_sym_table[i_seq];
00034 scitbx::vec3<double> const& site_i = sites_frac[i_seq];
00035 for(crystal::pair_sym_dict::const_iterator
00036 pair_sym_dict_i=pair_sym_dict.begin();
00037 pair_sym_dict_i!=pair_sym_dict.end();
00038 pair_sym_dict_i++) {
00039 unsigned j_seq = pair_sym_dict_i->first;
00040 scitbx::vec3<double> const& site_j = sites_frac[j_seq];
00041 af::const_ref<sgtbx::rt_mx> rt_mx_list = af::make_const_ref(
00042 pair_sym_dict_i->second);
00043 for(unsigned i_op=0;i_op<rt_mx_list.size();i_op++) {
00044 distances.push_back((
00045 orthogonalization_matrix
00046 * (site_i - rt_mx_list[i_op] * site_j)).length());
00047 }
00048 }
00049 }
00050 return distances;
00051 }
00052
00053 class adp_iso_local_sphere_restraints_energies
00054 {
00055 public:
00056 adp_iso_local_sphere_restraints_energies() {}
00057
00058 adp_iso_local_sphere_restraints_energies(
00059 af::const_ref<pair_sym_dict> const& pair_sym_table,
00060 scitbx::mat3<double> const& orthogonalization_matrix,
00061 af::const_ref<scitbx::vec3<double> > const& sites_frac,
00062 af::const_ref<double> const& u_isos,
00063 af::const_ref<bool> const& selection,
00064 af::const_ref<bool> const& use_u_iso,
00065 af::const_ref<bool> const& grad_u_iso,
00066 double sphere_radius,
00067 double distance_power,
00068 double average_power,
00069 double min_u_sum,
00070 bool compute_gradients,
00071 bool collect)
00072 {
00073 CCTBX_ASSERT(sites_frac.size() == pair_sym_table.size());
00074 CCTBX_ASSERT(u_isos.size() == pair_sym_table.size());
00075 CCTBX_ASSERT(use_u_iso.size() == pair_sym_table.size());
00076 CCTBX_ASSERT(grad_u_iso.size() == pair_sym_table.size());
00077 CCTBX_ASSERT(selection.size() == pair_sym_table.size());
00078 number_of_restraints = 0;
00079 residual_sum = 0;
00080 if (compute_gradients) {
00081 gradients.resize(u_isos.size());
00082 }
00083 for(unsigned i_seq=0;i_seq<pair_sym_table.size();i_seq++) {
00084 if(use_u_iso[i_seq] && selection[i_seq]) {
00085 crystal::pair_sym_dict const& pair_sym_dict = pair_sym_table[i_seq];
00086 scitbx::vec3<double> const& site_i = sites_frac[i_seq];
00087 for(crystal::pair_sym_dict::const_iterator
00088 pair_sym_dict_i=pair_sym_dict.begin();
00089 pair_sym_dict_i!=pair_sym_dict.end();
00090 pair_sym_dict_i++) {
00091 unsigned j_seq = pair_sym_dict_i->first;
00092 scitbx::vec3<double> const& site_j = sites_frac[j_seq];
00093 af::const_ref<sgtbx::rt_mx> rt_mx_list = af::make_const_ref(
00094 pair_sym_dict_i->second);
00095 for(unsigned i_op=0;i_op<rt_mx_list.size();i_op++) {
00096 double dist = (orthogonalization_matrix
00097 * (site_i - rt_mx_list[i_op] * site_j)).length();
00098 if (dist <= sphere_radius && dist > 0.0) {
00099 if(dist < 0.1) {
00100 dist = 0.1;
00101 }
00102 double one_over_weight = std::pow(dist, distance_power);
00103 CCTBX_ASSERT(one_over_weight != 0);
00104 double weight = 1./one_over_weight;
00105 double u1 = u_isos[i_seq];
00106 double u2 = u_isos[j_seq];
00107 if(u1 >= 0.0 && u2 >= 0) {
00108 double sum = u1 + u2;
00109 if (sum >= min_u_sum) {
00110 double ave_pow = std::pow(sum/2, average_power);
00111 CCTBX_ASSERT(ave_pow != 0);
00112 double u_diff = u1 - u2;
00113 double u_diff_sq = u_diff * u_diff;
00114 number_of_restraints++;
00115 residual_sum += (weight * u_diff_sq / ave_pow);
00116 if (compute_gradients && grad_u_iso[i_seq] &&
00117 grad_u_iso[j_seq]) {
00118 double mem1 = 2.* u_diff / ave_pow;
00119 CCTBX_ASSERT(sum * ave_pow != 0);
00120 double mem2 = average_power * u_diff_sq / (sum*ave_pow);
00121 gradients[i_seq] += weight * (mem1 - mem2);
00122 gradients[j_seq] += weight * (-mem1 - mem2);
00123 }
00124 if (collect) {
00125 u_i.push_back(u1);
00126 u_j.push_back(u2);
00127 r_ij.push_back(dist);
00128 }
00129 }
00130 }
00131 }
00132 }
00133 }
00134 }
00135 }
00136 }
00137
00138 unsigned number_of_restraints;
00139 double residual_sum;
00140 af::shared<double> gradients;
00141 af::shared<double> u_i;
00142 af::shared<double> u_j;
00143 af::shared<double> r_ij;
00144 };
00145
00153 inline
00154 af::shared<double>
00155 get_distances(
00156 af::const_ref<pair_sym_dict> const& pair_sym_table,
00157 af::const_ref<scitbx::vec3<double> > const& sites_cart)
00158 {
00159 CCTBX_ASSERT(sites_cart.size() == pair_sym_table.size());
00160 af::shared<double> distances;
00161 for(unsigned i_seq=0;i_seq<pair_sym_table.size();i_seq++) {
00162 crystal::pair_sym_dict const& pair_sym_dict = pair_sym_table[i_seq];
00163 scitbx::vec3<double> const& site_i = sites_cart[i_seq];
00164 for(crystal::pair_sym_dict::const_iterator
00165 pair_sym_dict_i=pair_sym_dict.begin();
00166 pair_sym_dict_i!=pair_sym_dict.end();
00167 pair_sym_dict_i++) {
00168 unsigned j_seq = pair_sym_dict_i->first;
00169 scitbx::vec3<double> const& site_j = sites_cart[j_seq];
00170 af::const_ref<sgtbx::rt_mx> rt_mx_list = af::make_const_ref(
00171 pair_sym_dict_i->second);
00172 CCTBX_ASSERT(rt_mx_list.size() == 1);
00173 CCTBX_ASSERT(rt_mx_list[0].is_unit_mx());
00174 distances.push_back((site_i - site_j).length());
00175 }
00176 }
00177 return distances;
00178 }
00179
00181 typedef std::set<unsigned> pair_asu_j_sym_group;
00183 typedef std::vector<pair_asu_j_sym_group> pair_asu_j_sym_groups;
00185 typedef std::map<unsigned, pair_asu_j_sym_groups> pair_asu_dict;
00187 typedef af::shared<pair_asu_dict> pair_asu_table_table;
00188
00192 template <typename FloatType=double, typename IntShiftType=int>
00193 class pair_asu_table
00194 {
00195 public:
00197 pair_asu_table() {}
00198
00200 pair_asu_table(
00201 boost::shared_ptr<
00202 direct_space_asu::asu_mappings<
00203 FloatType, IntShiftType> > asu_mappings)
00204 :
00205 asu_mappings_owner_(asu_mappings),
00206 asu_mappings_(asu_mappings.get()),
00207 table_(asu_mappings->mappings_const_ref().size())
00208 {}
00209
00211 void
00212 reserve(std::size_t n_sites_final) { table_.reserve(n_sites_final); }
00213
00215 void
00216 grow(std::size_t number_of_additional_entries)
00217 {
00218 table_.resize(table_.size()+number_of_additional_entries);
00219 }
00220
00222 boost::shared_ptr<
00223 direct_space_asu::asu_mappings<FloatType, IntShiftType> >
00224 asu_mappings() const { return asu_mappings_owner_; }
00225
00227 pair_asu_table_table const&
00228 table() const { return table_; }
00229
00231 bool
00232 contains(direct_space_asu::asu_mapping_index_pair const& pair) const
00233 {
00234 return contains(pair.i_seq, pair.j_seq, pair.j_sym);
00235 }
00236
00238 bool
00239 contains(unsigned i_seq, unsigned j_seq, unsigned j_sym) const
00240 {
00241 pair_asu_dict const& asu_dict = table_[i_seq];
00242 pair_asu_dict::const_iterator asu_dict_i = asu_dict.find(j_seq);
00243 if (asu_dict_i != asu_dict.end()) {
00244 for (pair_asu_j_sym_groups::const_iterator
00245 j_sym_groups_i = asu_dict_i->second.begin();
00246 j_sym_groups_i != asu_dict_i->second.end();
00247 j_sym_groups_i++) {
00248 if (j_sym_groups_i->find(j_sym) != j_sym_groups_i->end()) {
00249 return true;
00250 }
00251 }
00252 }
00253 return false;
00254 }
00255
00257 bool
00258 operator==(pair_asu_table const& other) const
00259 {
00260 af::const_ref<pair_asu_dict> tab_a = table_.const_ref();
00261 af::const_ref<pair_asu_dict> tab_b = other.table_.const_ref();
00262 if (tab_a.size() != tab_b.size()) return false;
00263 for(unsigned i_seq=0;i_seq<tab_a.size();i_seq++) {
00264 pair_asu_dict const& dict_a = tab_a[i_seq];
00265 pair_asu_dict const& dict_b = tab_b[i_seq];
00266 if (dict_a.size() != dict_b.size()) return false;
00267 for(pair_asu_dict::const_iterator
00268 dict_a_i=dict_a.begin();
00269 dict_a_i!=dict_a.end();
00270 dict_a_i++) {
00271 pair_asu_dict::const_iterator
00272 dict_b_i = dict_b.find(dict_a_i->first);
00273 if (dict_b_i == dict_b.end()) return false;
00274 pair_asu_j_sym_groups const& jgs_a = dict_a_i->second;
00275 pair_asu_j_sym_groups const& jgs_b = dict_b_i->second;
00276 if (jgs_a.size() != jgs_b.size()) return false;
00277 std::vector<bool> is_matched(jgs_b.size(), false);
00278 for(unsigned ig_a=0;ig_a<jgs_a.size();ig_a++) {
00279 unsigned ig_b=0;
00280 for(;ig_b<jgs_b.size();ig_b++) {
00281 if (is_matched[ig_b]) continue;
00282 if (jgs_a[ig_a] == jgs_b[ig_b]) {
00283 is_matched[ig_b] = true;
00284 break;
00285 }
00286 }
00287 if (ig_b == jgs_b.size()) return false;
00288 }
00289 }
00290 }
00291 return true;
00292 }
00293
00295 bool
00296 operator!=(pair_asu_table const& other) const
00297 {
00298 return !((*this) == other);
00299 }
00300
00302 af::shared<std::size_t>
00303 pair_counts() const
00304 {
00305 af::const_ref<pair_asu_dict> tab = table_.const_ref();
00306 af::shared<std::size_t> result((af::reserve(tab.size())));
00307 for(unsigned i_seq=0;i_seq<tab.size();i_seq++) {
00308 std::size_t count = 0;
00309 pair_asu_dict const& dict = tab[i_seq];
00310 for(pair_asu_dict::const_iterator
00311 dict_i=dict.begin();
00312 dict_i!=dict.end();
00313 dict_i++) {
00314 pair_asu_j_sym_groups const& jgs = dict_i->second;
00315 for(unsigned ig=0;ig<jgs.size();ig++) {
00316 count += jgs[ig].size();
00317 }
00318 }
00319 result.push_back(count);
00320 }
00321 return result;
00322 }
00323
00325
00328 af::shared<std::size_t>
00329 cluster_pivot_selection(
00330 bool general_positions_only=false,
00331 std::size_t max_clusters=0,
00332 unsigned estimated_reduction_factor=4) const
00333 {
00334 af::const_ref<pair_asu_dict> tab = table_.const_ref();
00335 af::shared<std::size_t> result;
00336 if (max_clusters > 0) {
00337 result.reserve(max_clusters);
00338 }
00339 else if (estimated_reduction_factor > 1) {
00340 result.reserve( (tab.size()+estimated_reduction_factor-1)
00341 / estimated_reduction_factor);
00342 }
00343 else {
00344 result.reserve(tab.size());
00345 }
00346 boost::scoped_array<bool> keep_flags(new bool[tab.size()]);
00347 for(unsigned i_seq=0;i_seq<tab.size();i_seq++) {
00348 if (general_positions_only
00349 && asu_mappings_->site_symmetry_table()
00350 .is_special_position(i_seq)) {
00351 keep_flags[i_seq] = false;
00352 continue;
00353 }
00354 bool keep = true;
00355 pair_asu_dict const& dict = tab[i_seq];
00356 for(pair_asu_dict::const_iterator
00357 dict_i=dict.begin();
00358 dict_i!=dict.end();
00359 dict_i++) {
00360 unsigned j_seq = dict_i->first;
00361 if (j_seq == i_seq || (j_seq < i_seq && keep_flags[j_seq])) {
00362 keep = false;
00363 break;
00364 }
00365 }
00366 keep_flags[i_seq] = keep;
00367 if (keep) {
00368 result.push_back(i_seq);
00369 if (result.size() == max_clusters) break;
00370 }
00371 }
00372 return result;
00373 }
00374
00380 pair_asu_table&
00381 add_all_pairs(
00382 FloatType const& distance_cutoff,
00383 FloatType const& epsilon=1.e-6)
00384 {
00385 neighbors::fast_pair_generator<FloatType, IntShiftType> pair_generator(
00386 asu_mappings_owner_,
00387 distance_cutoff*(1+epsilon),
00388 true);
00389 while (!pair_generator.at_end()) {
00390 add_pair(pair_generator.next());
00391 }
00392 return *this;
00393 }
00394
00396
00398 pair_asu_table&
00399 add_pair_sym_table(pair_sym_table const& sym_table)
00400 {
00401 af::const_ref<pair_sym_dict> sym_table_ref = sym_table.const_ref();
00402 for(unsigned i_seq=0;i_seq<sym_table_ref.size();i_seq++) {
00403 for(pair_sym_dict::const_iterator
00404 sym_dict_i = sym_table_ref[i_seq].begin();
00405 sym_dict_i != sym_table_ref[i_seq].end();
00406 sym_dict_i++) {
00407 unsigned j_seq = sym_dict_i->first;
00408 std::vector<sgtbx::rt_mx> const& rt_mx_list = sym_dict_i->second;
00409 for(unsigned i=0;i<rt_mx_list.size();i++) {
00410 add_pair(i_seq, j_seq, rt_mx_list[i]);
00411 }
00412 }
00413 }
00414 return *this;
00415 }
00416
00418 pair_asu_table&
00419 add_pair(direct_space_asu::asu_mapping_index_pair const& pair)
00420 {
00421 sgtbx::rt_mx rt_mx_i = asu_mappings_->get_rt_mx_i(pair);
00422 sgtbx::rt_mx rt_mx_j = asu_mappings_->get_rt_mx_j(pair);
00423 add_pair(
00424 pair.i_seq,
00425 pair.j_seq,
00426 rt_mx_i.inverse().multiply(rt_mx_j));
00427 return *this;
00428 }
00429
00431
00437 pair_asu_table&
00438 add_pair(unsigned i_seq, unsigned j_seq, sgtbx::rt_mx const& rt_mx_ji)
00439 {
00440 bool is_new = process_pair(i_seq, j_seq, rt_mx_ji);
00441 if (is_new && i_seq != j_seq) {
00442 is_new = process_pair(j_seq, i_seq, rt_mx_ji.inverse_cancel());
00443 CCTBX_ASSERT(is_new);
00444 }
00445 return *this;
00446 }
00447
00451 pair_asu_table&
00452 add_pair(af::tiny<unsigned, 2> const& i_seqs)
00453 {
00454 sgtbx::rt_mx rt_mx_ji(1, 1);
00455 bool is_new = process_pair(i_seqs[0], i_seqs[1], rt_mx_ji);
00456 if (is_new && i_seqs[0] != i_seqs[1]) {
00457 is_new = process_pair(i_seqs[1], i_seqs[0], rt_mx_ji);
00458 CCTBX_ASSERT(is_new);
00459 }
00460 return *this;
00461 }
00462
00464
00471 pair_sym_table
00472 extract_pair_sym_table(bool skip_j_seq_less_than_i_seq=true) const
00473 {
00474 pair_sym_table sym_table(asu_mappings_->mappings_const_ref().size());
00475 af::const_ref<pair_asu_dict> table_ref = table_.const_ref();
00476 for(unsigned i_seq=0;i_seq<table_.size();i_seq++) {
00477 sgtbx::rt_mx
00478 rt_mx_i_inv = asu_mappings_->get_rt_mx(i_seq, 0).inverse();
00479 pair_asu_dict const& asu_dict = table_ref[i_seq];
00480 pair_sym_dict& sym_dict = sym_table[i_seq];
00481 for(pair_asu_dict::const_iterator
00482 asu_dict_i = asu_dict.begin();
00483 asu_dict_i != asu_dict.end();
00484 asu_dict_i++) {
00485 unsigned j_seq = asu_dict_i->first;
00486 if (skip_j_seq_less_than_i_seq && j_seq < i_seq) continue;
00487 std::vector<sgtbx::rt_mx>& rt_mx_list = sym_dict[j_seq];
00488 for(pair_asu_j_sym_groups::const_iterator
00489 j_sym_groups_i = asu_dict_i->second.begin();
00490 j_sym_groups_i != asu_dict_i->second.end();
00491 j_sym_groups_i++) {
00492 unsigned j_sym = *(j_sym_groups_i->begin());
00493 sgtbx::rt_mx rt_mx_j = asu_mappings_->get_rt_mx(j_seq, j_sym);
00494 rt_mx_list.push_back(rt_mx_i_inv.multiply(rt_mx_j));
00495 }
00496 }
00497 }
00498 return sym_table;
00499 }
00500
00509 bool
00510 process_pair(
00511 unsigned i_seq,
00512 unsigned j_seq,
00513 sgtbx::rt_mx const& rt_mx_ji)
00514 {
00515 sgtbx::rt_mx rt_mx_i = asu_mappings_->get_rt_mx(i_seq, 0);
00516 sgtbx::rt_mx rt_mx_j = rt_mx_i.multiply(rt_mx_ji);
00517 int j_sym = asu_mappings_->find_i_sym(j_seq, rt_mx_j);
00518 return process_pair(i_seq, j_seq, rt_mx_ji, rt_mx_i, j_sym);
00519 }
00520
00529 bool
00530 process_pair(
00531 unsigned i_seq,
00532 unsigned j_seq,
00533 sgtbx::rt_mx const& rt_mx_ji,
00534 sgtbx::rt_mx const& rt_mx_i,
00535 int j_sym)
00536 {
00537 CCTBX_ASSERT(j_sym >= 0);
00538 if (contains(i_seq, j_seq, j_sym)) {
00539 return false;
00540 }
00541 table_[i_seq][j_seq].push_back(pair_asu_j_sym_group());
00542 pair_asu_j_sym_group& j_syms = table_[i_seq][j_seq].back();
00543 sgtbx::site_symmetry_table const&
00544 site_symmetry_table = asu_mappings_->site_symmetry_table();
00545 if ( i_seq != j_seq
00546 && !site_symmetry_table.is_special_position(i_seq)) {
00547 j_syms.insert(j_sym);
00548 }
00549 else {
00550 af::const_ref<sgtbx::rt_mx> const&
00551 site_symmetry_matrices = site_symmetry_table
00552 .get(i_seq).matrices().const_ref();
00553 for(unsigned i_mi=0;i_mi<site_symmetry_matrices.size();i_mi++) {
00554 sgtbx::rt_mx const& mi = site_symmetry_matrices[i_mi];
00555 if (i_seq == j_seq) {
00556 sgtbx::rt_mx rt_mx_j_eq
00557 = rt_mx_i.multiply(rt_mx_ji.multiply(mi).inverse_cancel());
00558 int j_sym_eq = asu_mappings_->find_i_sym(j_seq, rt_mx_j_eq);
00559 CCTBX_ASSERT(j_sym_eq >= 0);
00560 j_syms.insert(j_sym_eq);
00561 }
00562 sgtbx::rt_mx rt_mx_j_eq = rt_mx_i.multiply(mi.multiply(rt_mx_ji));
00563 int j_sym_eq = asu_mappings_->find_i_sym(j_seq, rt_mx_j_eq);
00564 CCTBX_ASSERT(j_sym_eq >= 0);
00565 j_syms.insert(j_sym_eq);
00566 }
00567 }
00568 return true;
00569 }
00570
00572 pair_asu_table
00573 angle_pair_asu_table() const
00574 {
00575 pair_asu_table result(asu_mappings_owner_);
00576 af::const_ref<pair_asu_dict> table_ref = table_.const_ref();
00577 typedef af::tiny<unsigned, 2> tu2;
00578 std::vector<tu2> pair_list;
00579 pair_list.reserve(16);
00580 for(unsigned i_seq=0;i_seq<table_.size();i_seq++) {
00581 pair_list.clear();
00582 pair_asu_dict const& asu_dict = table_ref[i_seq];
00583 for(pair_asu_dict::const_iterator
00584 asu_dict_i = asu_dict.begin();
00585 asu_dict_i != asu_dict.end();
00586 asu_dict_i++) {
00587 unsigned j_seq = asu_dict_i->first;
00588 for(pair_asu_j_sym_groups::const_iterator
00589 j_sym_groups_i = asu_dict_i->second.begin();
00590 j_sym_groups_i != asu_dict_i->second.end();
00591 j_sym_groups_i++) {
00592 for(pair_asu_j_sym_group::const_iterator
00593 j_sym_i = j_sym_groups_i->begin();
00594 j_sym_i != j_sym_groups_i->end();
00595 j_sym_i++) {
00596 pair_list.push_back(tu2(j_seq, *j_sym_i));
00597 }
00598 }
00599 }
00600 unsigned pair_list_size = static_cast<unsigned>(pair_list.size());
00601 for(unsigned i_jj1=0;i_jj1+1<pair_list_size;i_jj1++) {
00602 tu2 const& jj1 = pair_list[i_jj1];
00603 sgtbx::rt_mx rt_mx_jj1_inv = asu_mappings_->get_rt_mx(
00604 jj1[0], jj1[1]).inverse();
00605 for(unsigned i_jj2=i_jj1+1;i_jj2<pair_list_size;i_jj2++) {
00606 tu2 const& jj2 = pair_list[i_jj2];
00607 result.add_pair(
00608 jj1[0],
00609 jj2[0],
00610 rt_mx_jj1_inv.multiply(asu_mappings_->get_rt_mx(
00611 jj2[0], jj2[1])));
00612 }
00613 }
00614 }
00615 return result;
00616 }
00617
00618 protected:
00619 boost::shared_ptr<
00620 direct_space_asu::asu_mappings<FloatType, IntShiftType> >
00621 asu_mappings_owner_;
00622 const direct_space_asu::asu_mappings<FloatType, IntShiftType>*
00623 asu_mappings_;
00624 pair_asu_table_table table_;
00625 };
00626
00627 }}
00628
00629 #endif // CCTBX_CRYSTAL_PAIR_TABLES_H