xmap.h

//C Copyright (C) 2000-2006 Kevin Cowtan and University of York
//L
//L  This library is free software and is distributed under the terms
//L  and conditions of version 2.1 of the GNU Lesser General Public
//L  Licence (LGPL) with the following additional clause:
//L
//L     `You may also combine or link a "work that uses the Library" to
//L     produce a work containing portions of the Library, and distribute
//L     that work under terms of your choice, provided that you give
//L     prominent notice with each copy of the work that the specified
//L     version of the Library is used in it, and that you include or
//L     provide public access to the complete corresponding
//L     machine-readable source code for the Library including whatever
//L     changes were used in the work. (i.e. If you make changes to the
//L     Library you must distribute those, but you do not need to
//L     distribute source or object code to those portions of the work
//L     not covered by this licence.)'
//L
//L  Note that this clause grants an additional right and does not impose
//L  any additional restriction, and so does not affect compatibility
//L  with the GNU General Public Licence (GPL). If you wish to negotiate
//L  other terms, please contact the maintainer.
//L
//L  You can redistribute it and/or modify the library under the terms of
//L  the GNU Lesser General Public License as published by the Free Software
//L  Foundation; either version 2.1 of the License, or (at your option) any
//L  later version.
//L
//L  This library is distributed in the hope that it will be useful, but
//L  WITHOUT ANY WARRANTY; without even the implied warranty of
//L  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
//L  Lesser General Public License for more details.
//L
//L  You should have received a copy of the CCP4 licence and/or GNU
//L  Lesser General Public License along with this library; if not, write
//L  to the CCP4 Secretary, Daresbury Laboratory, Warrington WA4 4AD, UK.
//L  The GNU Lesser General Public can also be obtained by writing to the
//L  Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
//L  MA 02111-1307 USA


#ifndef CLIPPER_XMAP
#define CLIPPER_XMAP


#include "fftmap.h"
#include "fftmap_sparse.h"
#include "derivs.h"


namespace clipper
{

  class Xmap_cacheobj
  {
  public:
    class Key
    {
    public:
      Key( const Spgr_descr& spgr_descr, const Grid_sampling& grid ) :
        spgr_descr_(spgr_descr), grid_sampling_(grid) {}
      const Spgr_descr& spgr_descr() const { return spgr_descr_; }
      const Grid_sampling& grid_sampling() const { return grid_sampling_; }
    private:
      Spgr_descr spgr_descr_;
      Grid_sampling grid_sampling_;
    };

    Xmap_cacheobj( const Key& xmap_cachekey );  
    bool matches( const Key& xmap_cachekey ) const; 
    String format() const;  
    // data
    Key key;                         
    Grid_sampling xtl_grid;          
    Grid_range asu_grid;             
    Grid_range map_grid;             
    int nsym;                        // number of symops
    std::vector<unsigned char> asu;  
    std::vector<Isymop> isymop;      
    std::vector<int> du, dv, dw;     
    Array2d<unsigned char> symperm;  
    Mat33<> mat_grid_orth;           
    static Mutex mutex;              
  };



  class Xmap_base
  {
  public:
    enum FFTtype { Default, Normal, Sparse };  

    bool is_null() const;

    const Cell& cell() const { return cell_; }
    const Spacegroup& spacegroup() const { return spacegroup_; }
    const Grid_sampling& grid_sampling() const { return grid_sam_; }
    const Grid_range& grid_asu() const { return cacheref.data().asu_grid; }

    inline Coord_grid coord_of( const int& index ) const
      { return cacheref.data().map_grid.deindex( index ); }

    inline int index_of( const Coord_grid& coord ) const {
      if ( cacheref.data().asu_grid.in_grid( coord ) ) {
        const int i = cacheref.data().map_grid.index( coord );
        if ( asu[ i ] == 0 ) return i;
      }
      return -1;
    }
    Coord_grid to_map_unit( const Coord_grid& pos ) const
      { return pos.unit( grid_sam_ ); }

    const RTop<>& operator_orth_grid() const { return rt_orth_grid; }
    const RTop<>& operator_grid_orth() const { return rt_grid_orth; }

    inline Coord_orth coord_orth( const Coord_map& cm ) const
      { return Coord_orth( rt_grid_orth.rot()*cm ); }

    inline Coord_map coord_map( const Coord_orth& co ) const
      { return Coord_map ( rt_orth_grid.rot()*co ); }

    bool in_map( const Coord_grid& ) const { return true; }
    template<class I> bool in_map( const Coord_map& cm ) const { return true; }

    int multiplicity( const Coord_grid& pos ) const;


    class Map_reference_base
    {
    public:
      inline const Xmap_base& base_xmap() const { return *map_; }
      inline const int& index() const { return index_; }
      bool last() const { return ( index_ >= map_->map_grid.size() ); }
    protected:
      const Xmap_base* map_;
      int index_;
    };


    class Map_reference_index : public Map_reference_base
    {
    public:
      Map_reference_index() {}
      explicit Map_reference_index( const Xmap_base& map )
        { map_ = &map; index_=0; next(); }
      Map_reference_index( const Xmap_base& map, const Coord_grid& pos ) { map_ = &map; int sym; map_->find_sym( pos, index_, sym ); }
      inline Coord_grid coord() const
        { return map_->map_grid.deindex(index_); }
      inline const Coord_orth coord_orth() const
        { return Coord_orth( map_->rt_grid_orth.rot() * coord().coord_map() ); }
      inline Map_reference_index& set_coord( const Coord_grid& pos )
        { int sym; map_->find_sym( pos, index_, sym ); return *this; }
      inline Map_reference_index& next() {
        do {
          index_++; if ( last() ) break;
        } while ( map_->asu[index_] != 0 );
        return *this;
      }
      /* Use for e.g. peak search. Valid for -1 <= du/dv/dw <= 1 only.
         \param du/dv/dw Coordinate offset. \return Map index. */
      inline int index_offset(const int& du,const int& dv,const int& dw) const {
        int i = index_ + du*map_->du[0] + dv*map_->dv[0] + dw*map_->dw[0];
        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] ) ) ); }
        return i;
      }
      // inherited functions listed for documentation purposes
      //-- const Xmap_base& base_xmap() const;
      //-- const int& index() const;
      //-- bool last() const;
    };


    class Map_reference_coord : public Map_reference_base
    {
    public:
      Map_reference_coord() {}
      explicit Map_reference_coord( const Xmap_base& map )
        { map_ = &map; index_ = 0; next(); }
      Map_reference_coord( const Xmap_base& map, const Coord_grid& pos ) {
        map_ = &map;
        pos_ = pos;
        map_->find_sym( pos_, index_, sym_ );
      }
      inline const Coord_grid& coord() const { return pos_; }
      inline const Coord_orth coord_orth() const
        { return Coord_orth( map_->rt_grid_orth.rot() * coord().coord_map() ); }
      inline const int& sym() const { return sym_; }
      Map_reference_coord& set_coord( const Coord_grid& pos );

      inline Map_reference_coord& next() {
        sym_ = 0;
        do {
          index_++; if ( last() ) break;
        } while ( map_->asu[index_] != 0 );
        pos_ = map_->map_grid.deindex(index_);
        return *this;
      }
      // Increment u,v,w
      inline Map_reference_coord& next_u() { pos_.u()++; index_ += map_->du[sym_]; if (map_->asu[index_] != 0) edge(); return *this; }  
      inline Map_reference_coord& next_v() { pos_.v()++; index_ += map_->dv[sym_]; if (map_->asu[index_] != 0) edge(); return *this; }  
      inline Map_reference_coord& next_w() { pos_.w()++; index_ += map_->dw[sym_]; if (map_->asu[index_] != 0) edge(); return *this; }  
      inline Map_reference_coord& prev_u() { pos_.u()--; index_ -= map_->du[sym_]; if (map_->asu[index_] != 0) edge(); return *this; }  
      inline Map_reference_coord& prev_v() { pos_.v()--; index_ -= map_->dv[sym_]; if (map_->asu[index_] != 0) edge(); return *this; }  
      inline Map_reference_coord& prev_w() { pos_.w()--; index_ -= map_->dw[sym_]; if (map_->asu[index_] != 0) edge(); return *this; }  

      inline Map_reference_coord& operator =( const Coord_grid& pos )
        { return set_coord( pos ); }
      // inherited functions listed for documentation purposes
      //-- const Xmap_base& base_xmap() const;
      //-- const int& index() const;
      //-- bool last() const;

    protected:
      int sym_;
      Coord_grid pos_;

      void edge();
    };

    Map_reference_index first() const { return Map_reference_index( *this ); }
    Map_reference_coord first_coord() const { return Map_reference_coord( *this ); }
    static FFTtype& default_type() { return default_type_; }
  protected:
    ObjectCache<Xmap_cacheobj>::Reference cacheref;  
    const unsigned char* asu;  
    const Isymop* isymop;      
    const int* du;             
    const int* dv;             
    const int* dw;             
    Grid_range asu_grid;       
    Grid_range map_grid;       
    int nsym;                  

    Cell cell_;                    
    Spacegroup spacegroup_;        
    Grid_sampling grid_sam_;       

    RTop<> rt_orth_grid;           
    RTop<> rt_grid_orth;           

    Xmap_base();
    void init( const Spacegroup& spacegroup, const Cell& cell, const Grid_sampling& grid_sam );
    inline void find_sym( const Coord_grid& base, int& index, int& sym ) const;
    void asu_error( const Coord_grid& pos ) const;

    static FFTtype default_type_;       

    friend class Xmap_base::Map_reference_base;
    friend class Xmap_base::Map_reference_index;
    friend class Xmap_base::Map_reference_coord;
  };





  template<class T> class Xmap : public Xmap_base
  {
  public:
    Xmap() {}
    Xmap( const Spacegroup& spacegroup, const Cell& cell, const Grid_sampling& grid_sam ) { init( spacegroup, cell, grid_sam ); }
    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() ); }

    inline const T& operator[] (const Xmap_base::Map_reference_index& ix) const
      { return list[ix.index()]; }
    inline T& operator[] (const Xmap_base::Map_reference_index& ix)
      { return list[ix.index()]; }
    
    inline const T& operator[] (const Xmap_base::Map_reference_coord& ix) const
      { return list[ix.index()]; }
    inline T& operator[] (const Xmap_base::Map_reference_coord& ix)
      { return list[ix.index()]; }
    
    const T& get_data( const Coord_grid& pos ) const;
    void set_data( const Coord_grid& pos, const T& val );
    inline const T& get_data( const int& index ) const;
    bool set_data( const int& index, const T& val  );

    template<class I> T interp( const Coord_frac& pos ) const;
    template<class I> void interp_grad( const Coord_frac& pos, T& val, Grad_frac<T>& grad ) const;
    template<class I> void interp_curv( const Coord_frac& pos, T& val, Grad_frac<T>& grad, Curv_frac<T>& curv ) const;
    template<class I> T interp( const Coord_map& pos ) const;
    template<class I> void interp_grad( const Coord_map& pos, T& val, Grad_map<T>& grad ) const;
    template<class I> void interp_curv( const Coord_map& pos, T& val, Grad_map<T>& grad, Curv_map<T>& curv ) const;

    template<class H> void fft_from( const H& fphidata, const FFTtype type = Default );
    template<class H> void fft_to  ( H& fphidata, const FFTtype type = Default ) const;

    // inherited functions listed for documentation purposes
    //-- const Cell& cell() const;
    //-- const Spacegroup& spacegroup() const;
    //-- const Grid_sampling& grid_sampling() const;
    //-- const Grid_range& grid_asu() const;
    //-- inline Coord_grid coord_of( const int& index ) const;
    //-- inline int index_of( const Coord_grid& coord ) const;
    //-- const int multiplicity( const Coord_grid& pos ) const;
    //-- const Coord_grid to_map_unit( const Coord_grid& pos ) const;
    //-- const Map_reference_index first() const;
    //-- const Map_reference_coord first_coord() const;

    const T& operator =( const T& value );
    const Xmap<T>& operator +=( const Xmap<T>& other );
    const Xmap<T>& operator -=( const Xmap<T>& other );

  private:
    std::vector<T> list;
  };



  // implementations

  void Xmap_base::find_sym( const Coord_grid& base, int& index, int& sym ) const
  {
    // first deal with first symop, and cache for highest performance
    Coord_grid rot = to_map_unit( base );
    if ( asu_grid.in_grid( rot ) ) {
      index = map_grid.index( rot );
      if ( asu[ index ] == 0 ) {
        sym = 0;
      } else {
        sym = asu[ index ] - 1;
        index = map_grid.index( to_map_unit( base.transform(isymop[sym]) ) );
      }
    } else {
      // now deal with other symops
      for ( sym = 1; sym < nsym; sym++ ) {
        rot = to_map_unit( base.transform( isymop[sym] ) );
        if ( asu_grid.in_grid( rot ) ) {
          index = map_grid.index( rot );
          if ( asu[ index ] == 0 ) return;
        }
      }
      index = 0;  // redundent, to avoid compiler warnings
      asu_error( base );
    }
    return;
  }


  template<class T> const T& Xmap<T>::get_data( const Coord_grid& pos ) const
  {
    int index, sym;
    find_sym( pos, index, sym );
    return list[ index ];
  }

  template<class T> void Xmap<T>::set_data( const Coord_grid& pos, const T& val )
  {
    int index, sym;
    find_sym( pos, index, sym );
    list[ index ] = val;
  }

  template<class T> const T& Xmap<T>::get_data( const int& index ) const
  { return list[index]; }

  template<class T> bool Xmap<T>::set_data( const int& index, const T& val  )
  {
    if ( index >= 0 && index < list.size() )
      if ( asu[index] == 0 ) {
        list[index] = val;
        return true;
      }
    return false;
  }

  template<class T> template<class I> T Xmap<T>::interp( const Coord_frac& pos ) const
  {
    T val;
    I::interp( *this, pos.coord_map( grid_sam_ ), val );
    return val;
  }


  template<class T> template<class I> void Xmap<T>::interp_grad( const Coord_frac& pos, T& val, Grad_frac<T>& grad ) const
  {
    Grad_map<T> g;
    I::interp_grad( *this, pos.coord_map( grid_sam_ ), val, g );
    grad = g.grad_frac( grid_sam_ );
  }


  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
  {
    Grad_map<T> g;
    Curv_map<T> c;
    I::interp_curv( *this, pos.coord_map( grid_sam_ ), val, g, c );
    grad = g.grad_frac( grid_sam_ );
    curv = c.curv_frac( grid_sam_ );
  }


  template<class T> template<class I> T Xmap<T>::interp( const Coord_map& pos ) const
  {
    T val;
    I::interp( *this, pos, val );
    return val;
  }


  template<class T> template<class I> void Xmap<T>::interp_grad( const Coord_map& pos, T& val, Grad_map<T>& grad ) const
    { I::interp_grad( *this, pos, val, grad ); }


  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
    { I::interp_curv( *this, pos, val, grad, curv ); }


  template<class T> template<class H> void Xmap<T>::fft_from( const H& fphidata, const FFTtype type )
  {
    if ( type == Sparse || ( type == Default && default_type() == Sparse ) ) {
      // make a sparse fftmap
      FFTmap_sparse_p1_hx fftmap( grid_sampling() );
      // copy from reflection data
      typename H::HKL_reference_index ih;
      ffttype f, phi0, phi1;
      int sym;
      for ( ih = fphidata.first_data(); !ih.last(); fphidata.next_data( ih ) ) {
        f = fphidata[ih].f();
        if ( f != 0.0 ) {
          phi0 = fphidata[ih].phi();
          const HKL& hkl = ih.hkl();
          fftmap.set_hkl( hkl,
                          std::complex<ffttype>( f*cos(phi0), f*sin(phi0) ) );
          for ( sym = 1; sym < spacegroup_.num_primops(); sym++ ) {
            phi1 = phi0 + hkl.sym_phase_shift( spacegroup_.symop(sym) );
            fftmap.set_hkl( hkl.transform( isymop[sym] ),
                            std::complex<ffttype>( f*cos(phi1), f*sin(phi1) ) );
          }
        }
      }
      // require output ASU coords
      for ( Map_reference_index ix = first(); !ix.last(); ix.next() )
        fftmap.require_real_data( ix.coord() );
      // do fft
      fftmap.fft_h_to_x(1.0/cell().volume());
      // fill map ASU
      for ( Map_reference_index ix = first(); !ix.last(); ix.next() )
        (*this)[ix] = fftmap.real_data( ix.coord() );
    } else {
      // make a normal fftmap
      FFTmap_p1 fftmap( grid_sampling() );
      // copy from reflection data
      typename H::HKL_reference_index ih;
      ffttype f, phi0, phi1;
      int sym;
      for ( ih = fphidata.first_data(); !ih.last(); fphidata.next_data( ih ) ) {
        f = fphidata[ih].f();
        if ( f != 0.0 ) {
          phi0 = fphidata[ih].phi();
          const HKL& hkl = ih.hkl();
          fftmap.set_hkl( hkl,
                          std::complex<ffttype>( f*cos(phi0), f*sin(phi0) ) );
          for ( sym = 1; sym < spacegroup_.num_primops(); sym++ ) {
            phi1 = phi0 + hkl.sym_phase_shift( spacegroup_.symop(sym) );
            fftmap.set_hkl( hkl.transform( isymop[sym] ),
                            std::complex<ffttype>( f*cos(phi1), f*sin(phi1) ) );
          }
        }
      }
      // do fft
      fftmap.fft_h_to_x(1.0/cell().volume());
      // fill map ASU
      for ( Map_reference_index ix = first(); !ix.last(); ix.next() )
        (*this)[ix] = fftmap.real_data( ix.coord() );
    }
  }


  template<class T> template<class H> void Xmap<T>::fft_to  ( H& fphidata, const FFTtype type ) const
  {
    if ( type == Sparse || ( type == Default && default_type() == Sparse ) ) {
      // make a sparse fftmap
      FFTmap_sparse_p1_xh fftmap( grid_sampling() );
      // copy from map data
      ffttype f;
      int sym;
      for ( Map_reference_index ix = first(); !ix.last(); ix.next() ) {
        f = (*this)[ix];
        if ( f != 0.0 ) {
          fftmap.real_data( ix.coord() ) = f;
          for ( sym = 1; sym < cacheref.data().nsym; sym++ )
            fftmap.real_data(
              ix.coord().transform( isymop[sym] ).unit( grid_sam_ ) ) = f;
        }
      }
      // require output ASU coords
      typename H::HKL_reference_index ih;
      for ( ih = fphidata.first(); !ih.last(); ih.next() )
        fftmap.require_hkl( ih.hkl() );
      // do fft
      fftmap.fft_x_to_h(cell().volume());
      // fill data ASU
      for ( ih = fphidata.first(); !ih.last(); ih.next() ) {
        std::complex<ffttype> c = fftmap.get_hkl( ih.hkl() );
        fphidata[ih].f() = std::abs(c);
        fphidata[ih].phi() = std::arg(c);
      }
    } else {
      // make a normal fftmap
      FFTmap_p1 fftmap( grid_sampling() );
      // copy from map data
      ffttype f;
      int sym;
      for ( Map_reference_index ix = first(); !ix.last(); ix.next() ) {
        f = (*this)[ix];
        if ( f != 0.0 ) {
          fftmap.real_data( ix.coord() ) = f;
          for ( sym = 1; sym < cacheref.data().nsym; sym++ )
            fftmap.real_data(
              ix.coord().transform( isymop[sym] ).unit( grid_sam_ ) ) = f;
        }
      }
      // do fft
      fftmap.fft_x_to_h(cell().volume());
      // fill data ASU
      typename H::HKL_reference_index ih;
      for ( ih = fphidata.first(); !ih.last(); ih.next() ) {
        std::complex<ffttype> c = fftmap.get_hkl( ih.hkl() );
        fphidata[ih].f() = std::abs(c);
        fphidata[ih].phi() = std::arg(c);
      }
    }
  }


  template<class T> const T& Xmap<T>::operator =( const T& value )
  {
    // copy value into map
    for ( Map_reference_index im = first(); !im.last(); im.next() )
      list[im.index()] = value;
    return value;
  }


  template<class T> const Xmap<T>& Xmap<T>::operator +=( const Xmap<T>& other )
  {
    if ( spacegroup().hash() != other.spacegroup().hash() ||
         grid_sampling() != other.grid_sampling() )
      Message::message( Message_fatal( "Xmap: map mismatch in +=" ) );
    for ( Map_reference_index im = first(); !im.last(); im.next() )
      list[im.index()] += other[im];
    return (*this);
  }

  template<class T> const Xmap<T>& Xmap<T>::operator -=( const Xmap<T>& other )
  {
    if ( spacegroup().hash() != other.spacegroup().hash() ||
         grid_sampling() != other.grid_sampling() )
      Message::message( Message_fatal( "Xmap: map mismatch in -=" ) );
    for ( Map_reference_index im = first(); !im.last(); im.next() )
      list[im.index()] -= other[im];
    return (*this);
  }


} // namespace clipper

#endif