00001 #ifndef CCTBX_CRYSTAL_SITE_CLUSTER_ANALYSIS_H
00002 #define CCTBX_CRYSTAL_SITE_CLUSTER_ANALYSIS_H
00003
00004 #include <cctbx/crystal/direct_space_asu.h>
00005 #include <scitbx/cubicles.h>
00006
00007 namespace cctbx { namespace crystal {
00008
00009 template <typename FloatType=double, typename IntShiftType=int>
00010 class site_cluster_analysis
00011 {
00012 public:
00013 FloatType min_cross_distance;
00014 FloatType min_self_distance;
00015 bool general_positions_only;
00016 FloatType min_distance_sym_equiv;
00017 bool assert_min_distance_sym_equiv;
00018 unsigned estimated_reduction_factor;
00019
00020 protected:
00021 FloatType min_cross_distance_sq_;
00022 FloatType min_self_distance_sq_;
00023 FloatType cubicle_epsilon_;
00024 typedef
00025 direct_space_asu::asu_mappings<FloatType, IntShiftType>
00026 asu_mappings_t;
00027 boost::shared_ptr<asu_mappings_t> asu_mappings_owner_;
00028 asu_mappings_t* asu_mappings_;
00029 typedef std::vector<direct_space_asu::asu_mapping_index>
00030 cubicle_content_t;
00031 scitbx::cubicles<cubicle_content_t, FloatType> cubicles_;
00032 std::vector<std::size_t> registry_new_;
00033
00034 public:
00036 site_cluster_analysis() {}
00037
00038 site_cluster_analysis(
00039 sgtbx::space_group const& space_group,
00040 direct_space_asu::float_asu<FloatType> const& asu,
00041 FloatType const& min_cross_distance_,
00042 FloatType const& min_self_distance_=-1,
00043 bool general_positions_only_=false,
00044 unsigned estimated_reduction_factor_=4,
00045 FloatType const& asu_mappings_buffer_thickness=-1,
00046 FloatType const& min_cubicle_edge=5,
00047 FloatType const& cubicle_epsilon=-1)
00048 :
00049 min_cross_distance(min_cross_distance_),
00050 min_self_distance(min_self_distance_ >= 0
00051 ? min_self_distance_
00052 : min_cross_distance),
00053 general_positions_only(general_positions_only_),
00054 min_distance_sym_equiv(0.5),
00055 assert_min_distance_sym_equiv(true),
00056 estimated_reduction_factor(estimated_reduction_factor_),
00057 min_cross_distance_sq_(min_cross_distance*min_cross_distance),
00058 min_self_distance_sq_(min_self_distance*min_self_distance),
00059 cubicle_epsilon_(cubicle_epsilon >= 0
00060 ? cubicle_epsilon
00061 : asu.is_inside_epsilon()),
00062 asu_mappings_owner_(new asu_mappings_t(
00063 space_group,
00064 asu,
00065 asu_mappings_buffer_thickness >= 0
00066 ? asu_mappings_buffer_thickness
00067 : std::max(min_cross_distance, min_self_distance))),
00068 asu_mappings_(asu_mappings_owner_.get()),
00069 cubicles_(
00070 asu_mappings_->asu_buffer().box_min( true),
00071 asu_mappings_->asu_buffer().box_span( true),
00072 std::max(
00073 std::max(min_cross_distance, min_self_distance),
00074 min_cubicle_edge),
00075 cubicle_epsilon_)
00076 {
00077 CCTBX_ASSERT(min_cross_distance > 0);
00078 CCTBX_ASSERT(min_self_distance >= 0);
00079 CCTBX_ASSERT(asu_mappings_->buffer_thickness()
00080 >= std::max(min_cross_distance, min_self_distance));
00081 }
00082
00083 boost::shared_ptr<
00084 direct_space_asu::asu_mappings<FloatType, IntShiftType> >
00085 asu_mappings() const { return asu_mappings_owner_; }
00086
00087 void
00088 insert_fixed_site_frac(
00089 fractional<FloatType> const& original_site,
00090 sgtbx::site_symmetry_ops const& site_symmetry_ops)
00091 {
00092 registry_new_.clear();
00093 direct_space_asu::asu_mapping_index mi;
00094 mi.i_seq = asu_mappings_->mappings().size();
00095 asu_mappings_->process(original_site, site_symmetry_ops);
00096 typename asu_mappings_t::array_of_mappings_for_one_site const&
00097 mappings_i = asu_mappings_->mappings_const_ref().back();
00098 registry_new_.reserve(mappings_i.size());
00099 for(mi.i_sym=0; mi.i_sym<mappings_i.size(); mi.i_sym++) {
00100 std::size_t i1d_cub = cubicles_.ref.accessor()(
00101 cubicles_.i_cubicle(mappings_i[mi.i_sym].mapped_site()));
00102 cubicles_.ref[i1d_cub].push_back(mi);
00103 }
00104 }
00105
00106 void
00107 insert_fixed_site_frac(
00108 fractional<FloatType> const& original_site)
00109 {
00110 insert_fixed_site_frac(
00111 original_site,
00112 sgtbx::site_symmetry(
00113 asu_mappings_->asu().unit_cell(),
00114 asu_mappings_->space_group(),
00115 original_site,
00116 min_distance_sym_equiv,
00117 assert_min_distance_sym_equiv));
00118 }
00119
00120 protected:
00121 void
00122 discard_last_core()
00123 {
00124 asu_mappings_->discard_last();
00125 for(std::size_t i=registry_new_.size();i!=0;) {
00126 cubicles_.ref[registry_new_[--i]].pop_back();
00127 }
00128 registry_new_.clear();
00129 }
00130
00131 public:
00132 void
00133 discard_last()
00134 {
00135 if (registry_new_.size() == 0) {
00136 throw std::runtime_error(
00137 "site_cluster_analysis::discard_last() failure."
00138 " Potential problems are:\n"
00139 " - discard_last() called twice\n"
00140 " - insert_fixed_site_frac() called previously\n"
00141 " - the previous process_*() call returned false");
00142 }
00143 discard_last_core();
00144 }
00145
00146 bool
00147 process_site_frac(
00148 fractional<FloatType> const& original_site,
00149 sgtbx::site_symmetry_ops const& site_symmetry_ops)
00150 {
00151 registry_new_.clear();
00152 if (general_positions_only
00153 && !site_symmetry_ops.is_point_group_1()) return false;
00154 direct_space_asu::asu_mapping_index mi;
00155 mi.i_seq = asu_mappings_->mappings().size();
00156 asu_mappings_->process(original_site, site_symmetry_ops);
00157 sgtbx::rt_mx
00158 rt_mx_i_inverse = asu_mappings_->get_rt_mx(mi.i_seq, 0).inverse();
00159 af::const_ref<
00160 typename asu_mappings_t::array_of_mappings_for_one_site> const&
00161 mappings = asu_mappings_->mappings_const_ref();
00162 typename asu_mappings_t::array_of_mappings_for_one_site const&
00163 mappings_i = mappings.back();
00164 scitbx::vec3<unsigned> n_cub = cubicles_.ref.accessor();
00165 registry_new_.reserve(mappings_i.size());
00166 for(mi.i_sym=0; mi.i_sym<mappings_i.size(); mi.i_sym++) {
00167 scitbx::vec3<unsigned>
00168 i_cub = cubicles_.i_cubicle(mappings_i[mi.i_sym].mapped_site());
00169 scitbx::vec3<unsigned> j_cub_min;
00170 scitbx::vec3<unsigned> j_cub_max;
00171 j_cub_min[0] = (i_cub[0] == 0 ? 0 : i_cub[0]-1);
00172 j_cub_max[0] = (i_cub[0] == n_cub[0]-1 ? i_cub[0] : i_cub[0]+1);
00173 j_cub_min[1] = (i_cub[1] == 0 ? 0 : i_cub[1]-1);
00174 j_cub_max[1] = (i_cub[1] == n_cub[1]-1 ? i_cub[1] : i_cub[1]+1);
00175 j_cub_min[2] = (i_cub[2] == 0 ? 0 : i_cub[2]-1);
00176 j_cub_max[2] = (i_cub[2] == n_cub[2]-1 ? i_cub[2] : i_cub[2]+1);
00177 scitbx::vec3<unsigned> j_cub;
00178 for(j_cub[0]=j_cub_min[0];j_cub[0]<=j_cub_max[0];j_cub[0]++) {
00179 for(j_cub[1]=j_cub_min[1];j_cub[1]<=j_cub_max[1];j_cub[1]++) {
00180 for(j_cub[2]=j_cub_min[2];j_cub[2]<=j_cub_max[2];j_cub[2]++) {
00181 const cubicle_content_t* cub_j = &cubicles_.ref(j_cub);
00182 typename cubicle_content_t::const_iterator cub_ji;
00183 if (mi.i_sym == 0) {
00184 for(cub_ji=cub_j->begin();
00185 cub_ji!=cub_j->end();
00186 cub_ji++) {
00187 scitbx::vec3<FloatType> diff_vec =
00188 mappings[cub_ji->i_seq][cub_ji->i_sym].mapped_site()
00189 - mappings_i[0].mapped_site();
00190 if (diff_vec.length_sq() > min_cross_distance_sq_) continue;
00191 discard_last_core();
00192 return false;
00193 }
00194 }
00195 else {
00196 for(cub_ji=cub_j->begin();
00197 cub_ji!=cub_j->end();
00198 cub_ji++) {
00199 if (cub_ji->i_sym != 0) continue;
00200 scitbx::vec3<FloatType> diff_vec =
00201 mappings_i[mi.i_sym].mapped_site()
00202 - mappings[cub_ji->i_seq][0].mapped_site();
00203 if (diff_vec.length_sq() >
00204 (cub_ji->i_seq == mi.i_seq
00205 ? min_self_distance_sq_
00206 : min_cross_distance_sq_)) continue;
00207 discard_last_core();
00208 return false;
00209 }
00210 }
00211 }}}
00212 std::size_t i1d_cub = cubicles_.ref.accessor()(i_cub);
00213 cubicles_.ref[i1d_cub].push_back(mi);
00214 registry_new_.push_back(i1d_cub);
00215 }
00216 return true;
00217 }
00218
00219 bool
00220 process_site_frac(
00221 fractional<FloatType> const& original_site)
00222 {
00223 return process_site_frac(
00224 original_site,
00225 sgtbx::site_symmetry(
00226 asu_mappings_->asu().unit_cell(),
00227 asu_mappings_->space_group(),
00228 original_site,
00229 min_distance_sym_equiv,
00230 assert_min_distance_sym_equiv));
00231 }
00232
00234 std::size_t
00235 reserve_additional(std::size_t n, bool apply_estimated_reduction_factor)
00236 {
00237 if (apply_estimated_reduction_factor
00238 && estimated_reduction_factor > 1) {
00239 n = (n+estimated_reduction_factor-1)/estimated_reduction_factor;
00240 }
00241 asu_mappings_->reserve(n + asu_mappings_->mappings_const_ref().size());
00242 return n;
00243 }
00244
00245 af::shared<std::size_t>
00246 process_sites_frac(
00247 af::const_ref<scitbx::vec3<FloatType> > const& original_sites,
00248 sgtbx::site_symmetry_table const& site_symmetry_table,
00249 std::size_t max_clusters=0)
00250 {
00251 CCTBX_ASSERT(site_symmetry_table.indices_const_ref().size()
00252 == original_sites.size());
00253 #define CCTBX_CRYSTAL_SITE_CLUSTER_ANALYSIS_CONSTRUCT_SELECTION \
00254 af::shared<std::size_t> result; \
00255 if (max_clusters == 0) { \
00256 result.reserve(reserve_additional(original_sites.size(), true)); \
00257 } \
00258 else { \
00259 result.reserve(reserve_additional(max_clusters, false)); \
00260 }
00261 CCTBX_CRYSTAL_SITE_CLUSTER_ANALYSIS_CONSTRUCT_SELECTION
00262 for(std::size_t i=0;i<original_sites.size();i++) {
00263 if (process_site_frac(
00264 original_sites[i], site_symmetry_table.get(i))) {
00265 result.push_back(i);
00266 if (result.size() == max_clusters) break;
00267 }
00268 }
00269 return result;
00270 }
00271
00272 af::shared<std::size_t>
00273 process_sites_frac(
00274 af::const_ref<scitbx::vec3<FloatType> > const& original_sites,
00275 std::size_t max_clusters=0)
00276 {
00277 CCTBX_CRYSTAL_SITE_CLUSTER_ANALYSIS_CONSTRUCT_SELECTION
00278 for(std::size_t i=0;i<original_sites.size();i++) {
00279 if (process_site_frac(original_sites[i])) {
00280 result.push_back(i);
00281 if (result.size() == max_clusters) break;
00282 }
00283 }
00284 return result;
00285 }
00286
00287 af::shared<std::size_t>
00288 process_sites_cart(
00289 af::const_ref<scitbx::vec3<FloatType> > const& original_sites,
00290 sgtbx::site_symmetry_table const& site_symmetry_table,
00291 std::size_t max_clusters=0)
00292 {
00293 CCTBX_ASSERT(site_symmetry_table.indices_const_ref().size()
00294 == original_sites.size());
00295 CCTBX_CRYSTAL_SITE_CLUSTER_ANALYSIS_CONSTRUCT_SELECTION
00296 uctbx::unit_cell const& uc = asu_mappings_->asu().unit_cell();
00297 for(std::size_t i=0;i<original_sites.size();i++) {
00298 if (process_site_frac(
00299 uc.fractionalize(original_sites[i]),
00300 site_symmetry_table.get(i))) {
00301 result.push_back(i);
00302 if (result.size() == max_clusters) break;
00303 }
00304 }
00305 return result;
00306 }
00307
00308 af::shared<std::size_t>
00309 process_sites_cart(
00310 af::const_ref<scitbx::vec3<FloatType> > const& original_sites,
00311 std::size_t max_clusters=0)
00312 {
00313 CCTBX_CRYSTAL_SITE_CLUSTER_ANALYSIS_CONSTRUCT_SELECTION
00314 uctbx::unit_cell const& uc = asu_mappings_->asu().unit_cell();
00315 for(std::size_t i=0;i<original_sites.size();i++) {
00316 if (process_site_frac(uc.fractionalize(original_sites[i]))) {
00317 result.push_back(i);
00318 if (result.size() == max_clusters) break;
00319 }
00320 }
00321 return result;
00322 }
00323 };
00324
00325 }}
00326
00327 #endif // CCTBX_SITE_CLUSTER_ANALYSIS_H
