00001
00002 #ifndef CCTBX_MILLER_LOOKUP_UTILS_H
00003 #define CCTBX_MILLER_LOOKUP_UTILS_H
00004
00005 #include <scitbx/constants.h>
00006 #include <cmath>
00007 #include <cstdio>
00008 #include <iostream>
00009
00010 #include <cctbx/miller/sym_equiv.h>
00011 #include <cctbx/miller/match_indices.h>
00012 #include <cctbx/miller/asu.h>
00013 #include <cctbx/sgtbx/reciprocal_space_asu.h>
00014 #include <cctbx/sgtbx/space_group_type.h>
00015
00016 #include <map>
00017 #include <vector>
00018
00019 namespace cctbx { namespace miller{
00020 namespace lookup_utils{
00021
00022
00023 template <typename FloatType=double>
00024 class lookup_tensor
00025 {
00026 typedef std::map<cctbx::miller::index<>,
00027 std::size_t,
00028 cctbx::miller::fast_less_than<> > lookup_map_type;
00029 public:
00031 lookup_tensor() {}
00033 lookup_tensor(
00034 scitbx::af::const_ref< cctbx::miller::index<> > const& hkl,
00035 cctbx::sgtbx::space_group const& space_group,
00036 bool const& anomalous_flag)
00037 :
00038 space_group_(space_group),
00039 sgtype_(space_group_),
00040 asu_choice_(sgtype_),
00041 hkl_lookup_(),
00042 n_duplicates_( 0 ),
00043 anomalous_flag_(anomalous_flag),
00044 n_indices_( hkl.size() )
00045 {
00046
00047 for (unsigned ii=0;ii<hkl.size();ii++){
00048 cctbx::miller::asym_index asumap(space_group_,
00049 asu_choice_,
00050 hkl[ii]);
00051 cctbx::miller::index_table_layout_adaptor asu_target_hkl;
00052 asu_target_hkl = asumap.one_column(anomalous_flag_);
00053
00054 lookup_map_type::const_iterator l = hkl_lookup_.find( asu_target_hkl.h() );
00055 if (l == hkl_lookup_.end()) {
00056 hkl_lookup_[ asu_target_hkl.h() ] = ii;
00057 } else {
00058 n_duplicates_++;
00059 }
00060
00061
00062
00063 }
00064 }
00065
00066 long
00067 n_duplicates()
00068 {
00069 return( n_duplicates_ );
00070 }
00071
00072 long
00073 find_hkl( cctbx::miller::index<> const& target_hkl)
00074 {
00075 long hkl_location;
00076
00077 cctbx::miller::asym_index asumap(space_group_,
00078 asu_choice_,
00079 target_hkl);
00080 cctbx::miller::index_table_layout_adaptor asu_target_hkl;
00081 asu_target_hkl = asumap.one_column(anomalous_flag_);
00082 lookup_map_type::const_iterator
00083 l = hkl_lookup_.find( asu_target_hkl.h() );
00084 if (l == hkl_lookup_.end()) {
00085 hkl_location = -1;
00086 }
00087 else {
00088 hkl_location = l->second;
00089 }
00090 if (hkl_location >= n_indices_){
00091 hkl_location = -1;
00092 }
00093 return (hkl_location);
00094 }
00095
00096 scitbx::af::shared<long>
00097 find_hkl( scitbx::af::const_ref<cctbx::miller::index<> >
00098 const& target_hkl )
00099 {
00100 scitbx::af::shared< long > permutation_table( target_hkl.size(), -1 );
00101 for (unsigned ii=0;ii<target_hkl.size();ii++){
00102 permutation_table[ii] = find_hkl(target_hkl[ii]);
00103 }
00104 return(permutation_table);
00105 }
00106
00107 protected:
00108 int n_duplicates_;
00109 int n_indices_;
00110 cctbx::sgtbx::space_group space_group_;
00111 cctbx::sgtbx::space_group_type sgtype_;
00112 cctbx::sgtbx::reciprocal_space::asu asu_choice_;
00113 lookup_map_type hkl_lookup_;
00114 bool anomalous_flag_;
00115
00116 };
00117
00118
00119
00120
00121
00122
00123
00124 template <typename FloatType=double>
00125 class local_neighbourhood
00126 {
00127 public:
00129 local_neighbourhood(){}
00131 local_neighbourhood(
00132 scitbx::af::const_ref< cctbx::miller::index<> > const& hkl,
00133 cctbx::sgtbx::space_group const& space_group,
00134 bool const& anomalous_flag,
00135 long const& radius
00136 )
00137 :
00138 lookup_tensor_(hkl, space_group, anomalous_flag),
00139 radius_(radius)
00140 {
00141 SCITBX_ASSERT( hkl.size() > 0 );
00142 for (unsigned ii=0;ii<hkl.size();ii++){
00143 hkl_.push_back( hkl[ii] );
00144 }
00145 }
00146
00148 std::vector< unsigned >
00149 construct_neighbourhood( cctbx::miller::index<> const& center_hkl )
00150 {
00151 std::vector<unsigned> neighbourhood_lookup_tensor;
00152 long manhattan_distance;
00153 long ht,kt,lt,index;
00154 for (int hh = -radius_ ; hh< radius_ +1; hh++){
00155 for (int kk = -radius_ ; kk< radius_ +1; kk++){
00156 for (int ll = -radius_ ; ll< radius_ +1; ll++){
00157
00158 manhattan_distance = std::abs(hh) +
00159 std::abs(kk) +
00160 std::abs(ll);
00161 if (manhattan_distance<=radius_){
00162 if (manhattan_distance>0){
00163 ht = center_hkl[0] + hh;
00164 kt = center_hkl[1] + kk;
00165 lt = center_hkl[2] + ll;
00166 cctbx::miller::index<> picked_neighbour(ht,kt,lt);
00167 index = lookup_tensor_.find_hkl(picked_neighbour);
00168 if (index >=0){
00169 neighbourhood_lookup_tensor.push_back( unsigned(index) );
00170 }
00171 }
00172 }
00173
00174
00175 }
00176 }
00177 }
00178
00179 return(neighbourhood_lookup_tensor);
00180 }
00181
00182
00183 std::vector< unsigned >
00184 construct_neighbourhood( unsigned const& center_hkl )
00185 {
00186 SCITBX_ASSERT( hkl_.size() > center_hkl);
00187 return( construct_neighbourhood(hkl_[center_hkl]) );
00188 }
00189
00190
00191
00192
00193 scitbx::af::shared< std::vector< unsigned > >
00194 construct_neighbourhood(
00195 scitbx::af::shared< cctbx::miller::index<> > const& hkl )
00196 {
00197 scitbx::af::shared< std::vector< unsigned > > result;
00198 for(unsigned ii=0;ii<hkl.size();ii++){
00199 result.push_back( construct_neighbourhood( hkl[ii] ) );
00200 }
00201 return(result);
00202 }
00203
00204
00205 scitbx::af::shared< std::vector< unsigned > >
00206 construct_neighbourhood(
00207 scitbx::af::shared< long > const& hkl )
00208 {
00209 scitbx::af::shared< std::vector<unsigned> > result;
00210 for (unsigned ii=0;ii<hkl.size();ii++){
00211 if ( hkl[ii]>=0 ){
00212 result.push_back( construct_neighbourhood( hkl[ii] ) );
00213 }
00214 else{
00215 std::vector<unsigned> tmp;
00216 result.push_back( tmp );
00217 }
00218 }
00219 return( result );
00220
00221 }
00222
00223
00224
00225 scitbx::af::shared< std::vector< unsigned > >
00226 construct_neighbourhood()
00227 {
00228 scitbx::af::shared< std::vector<unsigned> > result;
00229 for (unsigned ii=0;ii<hkl_.size();ii++){
00230 result.push_back( construct_neighbourhood( hkl_[ii] ) );
00231 }
00232 return( result );
00233
00234 }
00235
00236 protected:
00237 lookup_tensor<FloatType> lookup_tensor_;
00238 scitbx::af::shared< cctbx::miller::index<> > hkl_;
00239
00240 long radius_;
00241
00242 };
00243
00244
00245
00246
00247
00248 template <typename FloatType=double>
00249 class neighbourhood_list
00250 {
00251 public:
00253 neighbourhood_list(){}
00255 neighbourhood_list(
00256 scitbx::af::const_ref< cctbx::miller::index<> > const& hkl,
00257 cctbx::sgtbx::space_group const& space_group,
00258 bool const& anomalous_flag,
00259 std::size_t const& radius)
00260 :
00261 generator_(hkl,space_group,anomalous_flag,radius)
00262 {
00263
00264 for (std::size_t ii=0;ii<hkl.size();ii++){
00265
00266 nb_list_.push_back(
00267 generator_.construct_neighbourhood( hkl[ii] )
00268 );
00269
00270 }
00271
00272 }
00273
00274 scitbx::af::shared< std::vector<unsigned> >
00275 neighbour_list(){ return(nb_list_); }
00276
00277 std::vector<unsigned>
00278 neighbour_list( unsigned const& this_hkl )
00279 {
00280 return( nb_list_[ this_hkl ] );
00281 }
00282
00283 protected:
00284 local_neighbourhood<FloatType> generator_;
00285 scitbx::af::shared< std::vector<unsigned> > nb_list_;
00286
00287 };
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298 template <typename FloatType = double >
00299 class local_area
00300 {
00301 public:
00302 local_area(){}
00303 local_area(
00304 scitbx::af::const_ref< cctbx::miller::index<> > const& hkl,
00305 scitbx::af::const_ref< bool > const& property,
00306 cctbx::sgtbx::space_group const& space_group,
00307 bool const& anomalous_flag,
00308 std::size_t const& radius,
00309 std::size_t const& depth,
00310 std::size_t const& at_least_this_number_of_neighbours)
00311 :
00312 depth_(depth),
00313 local_area_locator_( hkl, space_group, anomalous_flag, radius),
00314 used_( hkl.size(), 0 ),
00315 average_number_of_neighbours_( 0 )
00316 {
00317
00318 SCITBX_ASSERT( property.size() == hkl.size() );
00319
00320
00321 neighbour_list_ = local_area_locator_.construct_neighbourhood();
00322
00323
00324
00325 unsigned last_size=0, tmp_size;
00326 for (unsigned ii=0;ii<hkl.size();ii++){
00327
00328 bool done_searching=true;
00329 unsigned iter=0;
00330 unsigned count=0;
00331 std::vector<unsigned> tmp_vector;
00332 area_.push_back( tmp_vector );
00333
00334 if (property[ii]){
00335 tmp_vector.push_back(ii);
00336
00337 used_[ii]=1;
00338 done_searching=false;
00339 }
00340 last_size = tmp_vector.size();
00341
00342
00343
00344 while ( !done_searching ){
00345 tmp_size = tmp_vector.size();
00346
00347 for (unsigned item=last_size-1;item<tmp_size;item++){
00348 std::vector<unsigned> tmp_nbl;
00349 tmp_nbl = neighbour_list_[ tmp_vector[item] ];
00350 for (unsigned kk=0;kk<tmp_nbl.size();kk++){
00351
00352
00353 if (used_[ tmp_nbl[kk] ] == 0){
00354
00355
00356 tmp_vector.push_back( tmp_nbl[kk] );
00357 used_[ tmp_nbl[kk] ] = 1;
00358 if ( property[ tmp_nbl[kk] ] ){
00359 count++;
00360 }
00361 }
00362 }
00363 }
00364
00365 last_size = tmp_size;
00366
00367 if (iter>=depth-1){
00368 done_searching = true;
00369 }
00370 if (count>= at_least_this_number_of_neighbours){
00371 done_searching = true;
00372 }
00373
00374 iter++;
00375 }
00376
00377 average_number_of_neighbours_+=count;
00378
00379
00380 for ( unsigned jj=0;jj<tmp_vector.size();jj++){
00381
00382 used_[ tmp_vector[jj] ] = 0;
00383
00384 if (property[ tmp_vector[jj] ] ){
00385 area_[ii].push_back( tmp_vector[jj] );
00386 }
00387 }
00388
00389 }
00390 average_number_of_neighbours_/=FloatType( hkl.size() );
00391 }
00392
00393
00394
00395
00396 scitbx::af::shared< std::vector<unsigned> >
00397 area()
00398 {
00399 return( area_ );
00400 }
00401
00402 scitbx::af::shared< std::vector<unsigned> > area_;
00403
00404 protected:
00405 unsigned depth_;
00406 local_neighbourhood<FloatType> local_area_locator_;
00407 scitbx::af::shared< std::vector<unsigned> > neighbour_list_;
00408 scitbx::af::shared<int> used_;
00409 FloatType average_number_of_neighbours_;
00410
00411 };
00412
00413
00414
00415
00416
00417
00418
00419
00420 }}}
00421 #endif // CCTBX_MILLER_LOOKUP_UTILS_H