Clipper
xmap.h
1 
4 //C Copyright (C) 2000-2006 Kevin Cowtan and University of York
5 //L
6 //L This library is free software and is distributed under the terms
7 //L and conditions of version 2.1 of the GNU Lesser General Public
8 //L Licence (LGPL) with the following additional clause:
9 //L
10 //L `You may also combine or link a "work that uses the Library" to
11 //L produce a work containing portions of the Library, and distribute
12 //L that work under terms of your choice, provided that you give
13 //L prominent notice with each copy of the work that the specified
14 //L version of the Library is used in it, and that you include or
15 //L provide public access to the complete corresponding
16 //L machine-readable source code for the Library including whatever
17 //L changes were used in the work. (i.e. If you make changes to the
18 //L Library you must distribute those, but you do not need to
19 //L distribute source or object code to those portions of the work
20 //L not covered by this licence.)'
21 //L
22 //L Note that this clause grants an additional right and does not impose
23 //L any additional restriction, and so does not affect compatibility
24 //L with the GNU General Public Licence (GPL). If you wish to negotiate
25 //L other terms, please contact the maintainer.
26 //L
27 //L You can redistribute it and/or modify the library under the terms of
28 //L the GNU Lesser General Public License as published by the Free Software
29 //L Foundation; either version 2.1 of the License, or (at your option) any
30 //L later version.
31 //L
32 //L This library is distributed in the hope that it will be useful, but
33 //L WITHOUT ANY WARRANTY; without even the implied warranty of
34 //L MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
35 //L Lesser General Public License for more details.
36 //L
37 //L You should have received a copy of the CCP4 licence and/or GNU
38 //L Lesser General Public License along with this library; if not, write
39 //L to the CCP4 Secretary, Daresbury Laboratory, Warrington WA4 4AD, UK.
40 //L The GNU Lesser General Public can also be obtained by writing to the
41 //L Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
42 //L MA 02111-1307 USA
43 
44 
45 #ifndef CLIPPER_XMAP
46 #define CLIPPER_XMAP
47 
48 
49 #include "fftmap.h"
50 #include "fftmap_sparse.h"
51 #include "derivs.h"
52 
53 
54 namespace clipper
55 {
56 
58  {
59  public:
60  class Key
61  {
62  public:
63  Key( const Spgr_descr& spgr_descr, const Grid_sampling& grid ) :
64  spgr_descr_(spgr_descr), grid_sampling_(grid) {}
65  const Spgr_descr& spgr_descr() const { return spgr_descr_; }
66  const Grid_sampling& grid_sampling() const { return grid_sampling_; }
67  private:
68  Spgr_descr spgr_descr_;
69  Grid_sampling grid_sampling_;
70  };
71 
72  Xmap_cacheobj( const Key& xmap_cachekey );
73  bool matches( const Key& xmap_cachekey ) const;
74  String format() const;
75  // data
76  Key key;
80  int nsym; // number of symops
81  std::vector<unsigned char> asu;
82  std::vector<Isymop> isymop;
83  std::vector<int> du, dv, dw;
86  static Mutex mutex;
87  };
88 
89 
91 
100  class Xmap_base
101  {
102  public:
103  enum FFTtype { Default, Normal, Sparse };
104 
106  bool is_null() const;
107 
109  const Cell& cell() const { return cell_; }
111  const Spacegroup& spacegroup() const { return spacegroup_; }
113  const Grid_sampling& grid_sampling() const { return grid_sam_; }
115  const Grid_range& grid_asu() const { return cacheref.data().asu_grid; }
117 
118  inline Coord_grid coord_of( const int& index ) const
119  { return cacheref.data().map_grid.deindex( index ); }
121 
123  inline int index_of( const Coord_grid& coord ) const {
124  if ( cacheref.data().asu_grid.in_grid( coord ) ) {
125  const int i = cacheref.data().map_grid.index( coord );
126  if ( asu[ i ] == 0 ) return i;
127  }
128  return -1;
129  }
131  Coord_grid to_map_unit( const Coord_grid& pos ) const
132  { return pos.unit( grid_sam_ ); }
133 
135  const RTop<>& operator_orth_grid() const { return rt_orth_grid; }
137  const RTop<>& operator_grid_orth() const { return rt_grid_orth; }
139 
141  inline Coord_orth coord_orth( const Coord_map& cm ) const
142  { return Coord_orth( rt_grid_orth.rot()*cm ); }
144 
146  inline Coord_map coord_map( const Coord_orth& co ) const
147  { return Coord_map ( rt_orth_grid.rot()*co ); }
148 
150  bool in_map( const Coord_grid& ) const { return true; }
152  template<class I> bool in_map( const Coord_map& cm ) const { return true; }
153 
155  int multiplicity( const Coord_grid& pos ) const;
156 
158 
163  {
164  public:
166  inline const Xmap_base& base_xmap() const { return *map_; }
168  inline const int& index() const { return index_; }
170  bool last() const { return ( index_ >= map_->map_grid.size() ); }
171  protected:
173  const Xmap_base* map_;
175  int index_;
176  };
177 
179 
191  {
192  public:
196  explicit Map_reference_index( const Xmap_base& map )
197  { map_ = &map; index_=0; next(); }
199  Map_reference_index( const Xmap_base& map, const Coord_grid& pos ) { map_ = &map; int sym; map_->find_sym( pos, index_, sym ); }
201  inline Coord_grid coord() const
202  { return map_->map_grid.deindex(index_); }
204  inline const Coord_orth coord_orth() const
205  { return Coord_orth( map_->rt_grid_orth.rot() * coord().coord_map() ); }
208  { int sym; map_->find_sym( pos, index_, sym ); return *this; }
211  do {
212  index_++; if ( last() ) break;
213  } while ( map_->asu[index_] != 0 );
214  return *this;
215  }
217  /* Use for e.g. peak search. Valid for -1 <= du/dv/dw <= 1 only.
218  \param du/dv/dw Coordinate offset. \return Map index. */
219  inline int index_offset(const int& du,const int& dv,const int& dw) const {
220  int i = index_ + du*map_->du[0] + dv*map_->dv[0] + dw*map_->dw[0];
221  if ( map_->asu[i] != 0 ) { i = map_->map_grid.index( map_->to_map_unit( map_->map_grid.deindex(i).transform( map_->isymop[map_->asu[i]-1] ) ) ); }
222  return i;
223  }
224  // inherited functions listed for documentation purposes
225  //-- const Xmap_base& base_xmap() const;
226  //-- const int& index() const;
227  //-- bool last() const;
228  };
229 
231 
243  {
244  public:
248  explicit Map_reference_coord( const Xmap_base& map )
249  { map_ = &map; index_ = 0; next(); }
251  Map_reference_coord( const Xmap_base& map, const Coord_grid& pos ) {
252  map_ = &map;
253  pos_ = pos;
254  map_->find_sym( pos_, index_, sym_ );
255  }
257  inline const Coord_grid& coord() const { return pos_; }
259  inline const Coord_orth coord_orth() const
260  { return Coord_orth( map_->rt_grid_orth.rot() * coord().coord_map() ); }
262  inline const int& sym() const { return sym_; }
266 
268  sym_ = 0;
269  do {
270  index_++; if ( last() ) break;
271  } while ( map_->asu[index_] != 0 );
273  return *this;
274  }
275  // Increment u,v,w
276  inline Map_reference_coord& next_u() { pos_.u()++; index_ += map_->du[sym_]; if (map_->asu[index_] != 0) edge(); return *this; }
277  inline Map_reference_coord& next_v() { pos_.v()++; index_ += map_->dv[sym_]; if (map_->asu[index_] != 0) edge(); return *this; }
278  inline Map_reference_coord& next_w() { pos_.w()++; index_ += map_->dw[sym_]; if (map_->asu[index_] != 0) edge(); return *this; }
279  inline Map_reference_coord& prev_u() { pos_.u()--; index_ -= map_->du[sym_]; if (map_->asu[index_] != 0) edge(); return *this; }
280  inline Map_reference_coord& prev_v() { pos_.v()--; index_ -= map_->dv[sym_]; if (map_->asu[index_] != 0) edge(); return *this; }
281  inline Map_reference_coord& prev_w() { pos_.w()--; index_ -= map_->dw[sym_]; if (map_->asu[index_] != 0) edge(); return *this; }
282  inline Map_reference_coord& operator =( const Coord_grid& pos )
284  { return set_coord( pos ); }
285  // inherited functions listed for documentation purposes
286  //-- const Xmap_base& base_xmap() const;
287  //-- const int& index() const;
288  //-- bool last() const;
289 
290  protected:
292  int sym_;
295 
297  void edge();
298  };
299 
301  Map_reference_index first() const { return Map_reference_index( *this ); }
305  static FFTtype& default_type() { return default_type_; }
306  protected:
308  const unsigned char* asu;
309  const Isymop* isymop;
310  const int* du;
311  const int* dv;
312  const int* dw;
315  int nsym;
316 
320 
323 
325  Xmap_base();
327  void init( const Spacegroup& spacegroup, const Cell& cell, const Grid_sampling& grid_sam );
328  inline void find_sym( const Coord_grid& base, int& index, int& sym ) const;
329  void asu_error( const Coord_grid& pos ) const;
330 
331  static FFTtype default_type_;
332 
333  friend class Xmap_base::Map_reference_base;
334  friend class Xmap_base::Map_reference_index;
335  friend class Xmap_base::Map_reference_coord;
336  };
337 
338 
339 
340 
342 
356  template<class T> class Xmap : public Xmap_base
357  {
358  public:
360  Xmap() {}
362  Xmap( const Spacegroup& spacegroup, const Cell& cell, const Grid_sampling& grid_sam ) { init( spacegroup, cell, grid_sam ); }
364  void init( const Spacegroup& spacegroup, const Cell& cell, const Grid_sampling& grid_sam ) { Xmap_base::init( spacegroup, cell, grid_sam ); list.resize( cacheref.data().asu.size() ); }
365 
367  inline const T& operator[] (const Xmap_base::Map_reference_index& ix) const
368  { return list[ix.index()]; }
371  { return list[ix.index()]; }
372 
374  inline const T& operator[] (const Xmap_base::Map_reference_coord& ix) const
375  { return list[ix.index()]; }
378  { return list[ix.index()]; }
379 
381  const T& get_data( const Coord_grid& pos ) const;
383  void set_data( const Coord_grid& pos, const T& val );
385  inline const T& get_data( const int& index ) const;
387  bool set_data( const int& index, const T& val );
388 
390  template<class I> T interp( const Coord_frac& pos ) const;
392  template<class I> void interp_grad( const Coord_frac& pos, T& val, Grad_frac<T>& grad ) const;
394  template<class I> void interp_curv( const Coord_frac& pos, T& val, Grad_frac<T>& grad, Curv_frac<T>& curv ) const;
396  template<class I> T interp( const Coord_map& pos ) const;
398  template<class I> void interp_grad( const Coord_map& pos, T& val, Grad_map<T>& grad ) const;
400  template<class I> void interp_curv( const Coord_map& pos, T& val, Grad_map<T>& grad, Curv_map<T>& curv ) const;
401 
403  template<class H> void fft_from( const H& fphidata, const FFTtype type = Default );
405  template<class H> void fft_to ( H& fphidata, const FFTtype type = Default ) const;
406 
407  // inherited functions listed for documentation purposes
408  //-- const Cell& cell() const;
409  //-- const Spacegroup& spacegroup() const;
410  //-- const Grid_sampling& grid_sampling() const;
411  //-- const Grid_range& grid_asu() const;
412  //-- inline Coord_grid coord_of( const int& index ) const;
413  //-- inline int index_of( const Coord_grid& coord ) const;
414  //-- const int multiplicity( const Coord_grid& pos ) const;
415  //-- const Coord_grid to_map_unit( const Coord_grid& pos ) const;
416  //-- const Map_reference_index first() const;
417  //-- const Map_reference_coord first_coord() const;
418 
420  const T& operator =( const T& value );
422  const Xmap<T>& operator +=( const Xmap<T>& other );
424  const Xmap<T>& operator -=( const Xmap<T>& other );
425 
426  private:
427  std::vector<T> list;
428  };
429 
430 
431 
432  // implementations
433 
434  void Xmap_base::find_sym( const Coord_grid& base, int& index, int& sym ) const
435  {
436  // first deal with first symop, and cache for highest performance
437  Coord_grid rot = to_map_unit( base );
438  if ( asu_grid.in_grid( rot ) ) {
439  index = map_grid.index( rot );
440  if ( asu[ index ] == 0 ) {
441  sym = 0;
442  } else {
443  sym = asu[ index ] - 1;
444  index = map_grid.index( to_map_unit( base.transform(isymop[sym]) ) );
445  }
446  } else {
447  // now deal with other symops
448  for ( sym = 1; sym < nsym; sym++ ) {
449  rot = to_map_unit( base.transform( isymop[sym] ) );
450  if ( asu_grid.in_grid( rot ) ) {
451  index = map_grid.index( rot );
452  if ( asu[ index ] == 0 ) return;
453  }
454  }
455  index = 0; // redundent, to avoid compiler warnings
456  asu_error( base );
457  }
458  return;
459  }
460 
461 
466  template<class T> const T& Xmap<T>::get_data( const Coord_grid& pos ) const
467  {
468  int index, sym;
469  find_sym( pos, index, sym );
470  return list[ index ];
471  }
472 
477  template<class T> void Xmap<T>::set_data( const Coord_grid& pos, const T& val )
478  {
479  int index, sym;
480  find_sym( pos, index, sym );
481  list[ index ] = val;
482  }
483 
490  template<class T> const T& Xmap<T>::get_data( const int& index ) const
491  { return list[index]; }
492 
500  template<class T> bool Xmap<T>::set_data( const int& index, const T& val )
501  {
502  if ( index >= 0 && index < list.size() )
503  if ( asu[index] == 0 ) {
504  list[index] = val;
505  return true;
506  }
507  return false;
508  }
509 
518  template<class T> template<class I> T Xmap<T>::interp( const Coord_frac& pos ) const
519  {
520  T val;
521  I::interp( *this, pos.coord_map( grid_sam_ ), val );
522  return val;
523  }
524 
525 
533  template<class T> template<class I> void Xmap<T>::interp_grad( const Coord_frac& pos, T& val, Grad_frac<T>& grad ) const
534  {
535  Grad_map<T> g;
536  I::interp_grad( *this, pos.coord_map( grid_sam_ ), val, g );
537  grad = g.grad_frac( grid_sam_ );
538  }
539 
540 
550  template<class T> template<class I> void Xmap<T>::interp_curv( const Coord_frac& pos, T& val, Grad_frac<T>& grad, Curv_frac<T>& curv ) const
551  {
552  Grad_map<T> g;
553  Curv_map<T> c;
554  I::interp_curv( *this, pos.coord_map( grid_sam_ ), val, g, c );
555  grad = g.grad_frac( grid_sam_ );
556  curv = c.curv_frac( grid_sam_ );
557  }
558 
559 
568  template<class T> template<class I> T Xmap<T>::interp( const Coord_map& pos ) const
569  {
570  T val;
571  I::interp( *this, pos, val );
572  return val;
573  }
574 
575 
583  template<class T> template<class I> void Xmap<T>::interp_grad( const Coord_map& pos, T& val, Grad_map<T>& grad ) const
584  { I::interp_grad( *this, pos, val, grad ); }
585 
586 
596  template<class T> template<class I> void Xmap<T>::interp_curv( const Coord_map& pos, T& val, Grad_map<T>& grad, Curv_map<T>& curv ) const
597  { I::interp_curv( *this, pos, val, grad, curv ); }
598 
599 
604  template<class T> template<class H> void Xmap<T>::fft_from( const H& fphidata, const FFTtype type )
605  {
606  if ( type == Sparse || ( type == Default && default_type() == Sparse ) ) {
607  // make a sparse fftmap
608  FFTmap_sparse_p1_hx fftmap( grid_sampling() );
609  // copy from reflection data
610  typename H::HKL_reference_index ih;
611  ffttype f, phi0, phi1;
612  int sym;
613  for ( ih = fphidata.first_data(); !ih.last(); fphidata.next_data( ih ) ) {
614  f = fphidata[ih].f();
615  if ( f != 0.0 ) {
616  phi0 = fphidata[ih].phi();
617  const HKL& hkl = ih.hkl();
618  fftmap.set_hkl( hkl,
619  std::complex<ffttype>( f*cos(phi0), f*sin(phi0) ) );
620  for ( sym = 1; sym < spacegroup_.num_primops(); sym++ ) {
621  phi1 = phi0 + hkl.sym_phase_shift( spacegroup_.symop(sym) );
622  fftmap.set_hkl( hkl.transform( isymop[sym] ),
623  std::complex<ffttype>( f*cos(phi1), f*sin(phi1) ) );
624  }
625  }
626  }
627  // require output ASU coords
628  for ( Map_reference_index ix = first(); !ix.last(); ix.next() )
629  fftmap.require_real_data( ix.coord() );
630  // do fft
631  fftmap.fft_h_to_x(1.0/cell().volume());
632  // fill map ASU
633  for ( Map_reference_index ix = first(); !ix.last(); ix.next() )
634  (*this)[ix] = fftmap.real_data( ix.coord() );
635  } else {
636  // make a normal fftmap
637  FFTmap_p1 fftmap( grid_sampling() );
638  // copy from reflection data
639  typename H::HKL_reference_index ih;
640  ffttype f, phi0, phi1;
641  int sym;
642  for ( ih = fphidata.first_data(); !ih.last(); fphidata.next_data( ih ) ) {
643  f = fphidata[ih].f();
644  if ( f != 0.0 ) {
645  phi0 = fphidata[ih].phi();
646  const HKL& hkl = ih.hkl();
647  fftmap.set_hkl( hkl,
648  std::complex<ffttype>( f*cos(phi0), f*sin(phi0) ) );
649  for ( sym = 1; sym < spacegroup_.num_primops(); sym++ ) {
650  phi1 = phi0 + hkl.sym_phase_shift( spacegroup_.symop(sym) );
651  fftmap.set_hkl( hkl.transform( isymop[sym] ),
652  std::complex<ffttype>( f*cos(phi1), f*sin(phi1) ) );
653  }
654  }
655  }
656  // do fft
657  fftmap.fft_h_to_x(1.0/cell().volume());
658  // fill map ASU
659  for ( Map_reference_index ix = first(); !ix.last(); ix.next() )
660  (*this)[ix] = fftmap.real_data( ix.coord() );
661  }
662  }
663 
664 
673  template<class T> template<class H> void Xmap<T>::fft_to ( H& fphidata, const FFTtype type ) const
674  {
675  if ( type == Sparse || ( type == Default && default_type() == Sparse ) ) {
676  // make a sparse fftmap
677  FFTmap_sparse_p1_xh fftmap( grid_sampling() );
678  // copy from map data
679  ffttype f;
680  int sym;
681  for ( Map_reference_index ix = first(); !ix.last(); ix.next() ) {
682  f = (*this)[ix];
683  if ( f != 0.0 ) {
684  fftmap.real_data( ix.coord() ) = f;
685  for ( sym = 1; sym < cacheref.data().nsym; sym++ )
686  fftmap.real_data(
687  ix.coord().transform( isymop[sym] ).unit( grid_sam_ ) ) = f;
688  }
689  }
690  // require output ASU coords
691  typename H::HKL_reference_index ih;
692  for ( ih = fphidata.first(); !ih.last(); ih.next() )
693  fftmap.require_hkl( ih.hkl() );
694  // do fft
695  fftmap.fft_x_to_h(cell().volume());
696  // fill data ASU
697  for ( ih = fphidata.first(); !ih.last(); ih.next() ) {
698  std::complex<ffttype> c = fftmap.get_hkl( ih.hkl() );
699  fphidata[ih].f() = std::abs(c);
700  fphidata[ih].phi() = std::arg(c);
701  }
702  } else {
703  // make a normal fftmap
704  FFTmap_p1 fftmap( grid_sampling() );
705  // copy from map data
706  ffttype f;
707  int sym;
708  for ( Map_reference_index ix = first(); !ix.last(); ix.next() ) {
709  f = (*this)[ix];
710  if ( f != 0.0 ) {
711  fftmap.real_data( ix.coord() ) = f;
712  for ( sym = 1; sym < cacheref.data().nsym; sym++ )
713  fftmap.real_data(
714  ix.coord().transform( isymop[sym] ).unit( grid_sam_ ) ) = f;
715  }
716  }
717  // do fft
718  fftmap.fft_x_to_h(cell().volume());
719  // fill data ASU
720  typename H::HKL_reference_index ih;
721  for ( ih = fphidata.first(); !ih.last(); ih.next() ) {
722  std::complex<ffttype> c = fftmap.get_hkl( ih.hkl() );
723  fphidata[ih].f() = std::abs(c);
724  fphidata[ih].phi() = std::arg(c);
725  }
726  }
727  }
728 
729 
732  template<class T> const T& Xmap<T>::operator =( const T& value )
733  {
734  // copy value into map
735  for ( Map_reference_index im = first(); !im.last(); im.next() )
736  list[im.index()] = value;
737  return value;
738  }
739 
740 
742  template<class T> const Xmap<T>& Xmap<T>::operator +=( const Xmap<T>& other )
743  {
744  if ( spacegroup().hash() != other.spacegroup().hash() ||
745  grid_sampling() != other.grid_sampling() )
746  Message::message( Message_fatal( "Xmap: map mismatch in +=" ) );
747  for ( Map_reference_index im = first(); !im.last(); im.next() )
748  list[im.index()] += other[im];
749  return (*this);
750  }
751 
753  template<class T> const Xmap<T>& Xmap<T>::operator -=( const Xmap<T>& other )
754  {
755  if ( spacegroup().hash() != other.spacegroup().hash() ||
756  grid_sampling() != other.grid_sampling() )
757  Message::message( Message_fatal( "Xmap: map mismatch in -=" ) );
758  for ( Map_reference_index im = first(); !im.last(); im.next() )
759  list[im.index()] -= other[im];
760  return (*this);
761  }
762 
763 
764 } // namespace clipper
765 
766 #endif