hkl_data.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_HKL_DATA
#define CLIPPER_HKL_DATA


#include "hkl_info.h"


namespace clipper
{

  class HKL_data_cacheobj : public HKL_info
  {
  public:
    class Key
    {
    public:
      Key( const Spgr_descr& spgr_descr, const Cell& cell_descr, const HKL_sampling& hkl_sam ) : spgr_descr_(spgr_descr), cell_descr_(cell_descr), hkl_sampling_(hkl_sam) {}
      const Spgr_descr& spgr_descr() const { return spgr_descr_; }
      const Cell_descr& cell_descr() const { return cell_descr_; }
      const HKL_sampling& hkl_sampling() const { return hkl_sampling_; }
    private:
      Spgr_descr spgr_descr_;
      Cell_descr cell_descr_;
      HKL_sampling hkl_sampling_;
    };

    HKL_data_cacheobj( const Key& hkl_data_cachekey );
    bool matches( const Key& hkl_data_cachekey ) const;
    String format() const;
    static Mutex mutex;                 
  private:
    Key key;
  };



  class Datatype_base
  {
  protected:
    //-- Datatype_base();
    void set_null();
    static String type();
    void friedel();
    void shift_phase(const ftype& dphi);
    bool missing() const;
    static int data_size();
    static String data_names();
    void data_export( xtype array[] ) const;
    void data_import( const xtype array[] );
  };



  class HKL_data_base
  {
   public:
    // Coordinate reference types
    typedef HKL_info::HKL_reference_index HKL_reference_index;
    typedef HKL_info::HKL_reference_coord HKL_reference_coord;

    virtual void init( const HKL_info& hkl_info, const Cell& cell );
    virtual void init( const HKL_data_base& hkl_data );
    virtual void init( const Spacegroup& spacegroup, const Cell& cell, const HKL_sampling& hkl_sampling );

    // generic methods
    bool is_null() const;

    const HKL_info& base_hkl_info() const { return *parent_hkl_info; }
    const Cell& base_cell() const { return *parent_cell; }

    const Spacegroup& spacegroup() const { return spacegroup_; }
    const Cell& cell() const { return cell_; }
    const Resolution& resolution() const { return resolution_; }
    const HKL_sampling& hkl_sampling() const { return hkl_sampling_; }
    const HKL_info& hkl_info() const { return *parent_hkl_info; }

    ftype invresolsq( const int& index ) const;
    Range<ftype> invresolsq_range() const;
    int num_obs() const;

    virtual void update() = 0;
    virtual String type() const = 0;
    virtual bool missing(const int& index) const = 0;
    virtual void set_null(const int& index) = 0;
    virtual int data_size() const = 0;
    virtual String data_names() const = 0;
    virtual void data_export( const HKL& hkl, xtype array[] ) const = 0;
    virtual void data_import( const HKL& hkl, const xtype array[] ) = 0;
    virtual void mask(const HKL_data_base& mask) = 0;

    HKL_reference_index first() const;
    HKL_reference_index first_data() const;
    HKL_reference_index& next_data( HKL_reference_index& ih ) const;

    void debug() const;

   protected:
    const HKL_info* parent_hkl_info;
    const Cell* parent_cell;
    bool cell_matches_parent;

    // clipper2 members
    ObjectCache<HKL_data_cacheobj>::Reference cacheref;  
    Spacegroup spacegroup_;
    Cell cell_;
    HKL_sampling hkl_sampling_;
    Resolution resolution_;

    HKL_data_base();
    virtual ~HKL_data_base() {}
  };



  template<class T> class HKL_data : public HKL_data_base
  {
  public:
    HKL_data() {}
    explicit HKL_data( const HKL_info& hkl_info );
    HKL_data( const HKL_info& hkl_info, const Cell& cell );
    HKL_data( const Spacegroup& spacegroup, const Cell& cell, const HKL_sampling& hkl_sampling );
    explicit HKL_data( const HKL_data_base& hkl_data );

    void init( const HKL_info& hkl_info, const Cell& cell );
    void init( const Spacegroup& spacegroup, const Cell& cell, const HKL_sampling& hkl_sampling );
    void init( const HKL_data_base& hkl_data );
    void update();

    // type specific methods
    String type() const { return T::type(); }
    bool missing(const int& index) const { return list[index].missing(); }
    void set_null(const int& index) { list[index].set_null(); }
    int data_size() const { return T::data_size(); }
    String data_names() const { return T::data_names(); }
    void data_export( const HKL& hkl, xtype array[] ) const
      { T datum; get_data( hkl, datum ); datum.data_export( array ); }
    void data_import( const HKL& hkl, const xtype array[] )
      { T datum; datum.data_import( array ); set_data( hkl, datum ); }
    void mask(const HKL_data_base& mask);

    // data access methods: by HKL_reference_index
    const T& operator[] (const HKL_info::HKL_reference_index& i) const
      { return list[i.index()]; }
    T& operator[] (const HKL_info::HKL_reference_index& i)
      { return list[i.index()]; }

    // data access methods: by HKL_reference_coord
    T operator[] (const HKL_info::HKL_reference_coord& ih) const;
    bool get_data(const HKL_info::HKL_reference_coord& ih, T& data) const;
    bool set_data(const HKL_info::HKL_reference_coord& ih, const T& data);

    // data access methods: by index
    const T& operator[] (const int& index) const { return list[index]; }
    T& operator[] (const int& index) { return list[index]; }

    // data access methods: by hkl
    T operator[] (const HKL& hkl) const;
    bool get_data(const HKL& hkl, T& data) const;
    bool set_data(const HKL& hkl, const T& data);

    // COMPUTATION OPERATORS
    template<class C> void compute( const C& op )
      { for (HKL_info::HKL_reference_index ih=parent_hkl_info->first(); !ih.last(); ih.next()) list[ih.index()] = op( ih ); }
    template<class S, class C> void compute( const HKL_data<S>& src, const C& op )
      { for (HKL_info::HKL_reference_index ih=parent_hkl_info->first(); !ih.last(); ih.next()) list[ih.index()] = op( ih, src[ih] ); }
    template<class S1, class S2, class C> void compute( const HKL_data<S1>& src1, const HKL_data<S2>& src2, const C& op )
      { for (HKL_info::HKL_reference_index ih=parent_hkl_info->first(); !ih.last(); ih.next()) list[ih.index()] = op( ih, src1[ih], src2[ih] ); }

    // inherited functions lists for documentation purposes
    //-- const HKL_info& base_hkl_info() const;
    //-- const Cell& base_cell() const;
    //-- const ftype invresolsq(const int& index) const;
    //-- const Range<ftype> invresolsq_range() const;
    //-- const int num_obs() const;
    //-- HKL_reference_index first() const;
    //-- HKL_reference_index first_data() const;
    //-- HKL_reference_index& next_data( HKL_reference_index& ih ) const;

    HKL_data<T>& operator =( const HKL_data<T>& other );
    HKL_data<T>& operator =( const T& value );

    void debug() const;

  protected:
    // members
    std::vector<T> list;
  };







  // Template implementations

  ftype HKL_info::HKL_reference_base::invresolsq( const HKL_data_base& hkldata ) const
    { return hkldata.invresolsq( index_ ); }

  template<class T> HKL_data<T>::HKL_data( const HKL_info& hkl_info )
  {
    init( hkl_info, hkl_info.cell() );
  }

  template<class T> HKL_data<T>::HKL_data( const HKL_info& hkl_info, const Cell& cell )
  {
    init( hkl_info, cell );
  }

  template<class T> HKL_data<T>::HKL_data( const Spacegroup& spacegroup, const Cell& cell, const HKL_sampling& hkl_sampling )
  {
    init( spacegroup, cell, hkl_sampling );
  }

  template<class T> HKL_data<T>::HKL_data( const HKL_data_base& hkl_data )
  {
    init( hkl_data );
  }

  template<class T> void HKL_data<T>::init( const HKL_info& hkl_info, const Cell& cell )
  {
    HKL_data_base::init( hkl_info, cell );
    update();
  }

  template<class T> void HKL_data<T>::init( const Spacegroup& spacegroup, const Cell& cell, const HKL_sampling& hkl_sampling )
  {
    HKL_data_base::init( spacegroup, cell, hkl_sampling );
    update();
  }

  template<class T> void HKL_data<T>::init( const HKL_data_base& hkl_data )
  {
    HKL_data_base::init( hkl_data );
    update();
  }

  template<class T> void HKL_data<T>::update()
  {
    if ( parent_hkl_info != NULL ) {
      T null; null.set_null();
      list.resize( parent_hkl_info->num_reflections(), null );
    }
  }

  template<class T> void HKL_data<T>::mask(const HKL_data_base& mask)
  {
    T null; null.set_null();
    for ( int i = 0; i < list.size(); i++ )
      if ( mask.missing(i) ) list[i] = null;
  }

  template<class T> T HKL_data<T>::operator[] (const HKL_info::HKL_reference_coord& ih) const
  {
    if ( ih.index() < 0 ) { T null; null.set_null(); return null; }
    T data = list[ih.index()];
    if ( ih.friedel() ) data.friedel();
    data.shift_phase( -ih.hkl().sym_phase_shift( parent_hkl_info->spacegroup().symop(ih.sym()) ) );
    return data;
  }

  template<class T> bool HKL_data<T>::get_data(const HKL_info::HKL_reference_coord& ih, T& data) const
  {
    if ( ih.index() < 0 ) { data.set_null(); return false; }
    data = list[ih.index()];
    if ( ih.friedel() ) data.friedel();
    data.shift_phase( -ih.hkl().sym_phase_shift( parent_hkl_info->spacegroup().symop(ih.sym()) ) );
    return true;
  }

  template<class T> bool HKL_data<T>::set_data(const HKL_info::HKL_reference_coord& ih, const T& data)
  {
    if ( ih.index() < 0 ) return false;
    T& ldata = list[ih.index()];
    ldata = data;
    ldata.shift_phase( ih.hkl().sym_phase_shift( parent_hkl_info->spacegroup().symop(ih.sym()) ) );
    if ( ih.friedel() ) ldata.friedel();
    return true;
  }

  template<class T> T HKL_data<T>::operator[] (const HKL& hkl) const
  {
    int index, sym; bool friedel;

    index = parent_hkl_info->index_of( parent_hkl_info->
                                       find_sym(hkl, sym, friedel) );
    if ( index < 0 ) { T null; null.set_null(); return null; }
    T data = list[index];
    if (friedel) data.friedel();
    data.shift_phase( -hkl.sym_phase_shift( parent_hkl_info->spacegroup().symop(sym) ) );
    return data;
  }

  template<class T> bool HKL_data<T>::get_data(const HKL& hkl, T& data) const
  {
    int index, sym; bool friedel;

    index = parent_hkl_info->index_of( parent_hkl_info->
                                        find_sym(hkl, sym, friedel) );
    if ( index < 0 ) { data.set_null(); return false; }
    data = list[index];
    if (friedel) data.friedel();
    data.shift_phase( -hkl.sym_phase_shift(parent_hkl_info->spacegroup().symop(sym)) );
    return true;
  }

  template<class T> bool HKL_data<T>::set_data(const HKL& hkl, const T& data_)
  {
    int index, sym; bool friedel;
    index = parent_hkl_info->index_of( parent_hkl_info->
                                       find_sym(hkl, sym, friedel) );
    if ( index < 0 ) { return false; }
    T& ldata = list[index];
    ldata = data_;
    ldata.shift_phase( hkl.sym_phase_shift(parent_hkl_info->spacegroup().symop(sym)) );
    if (friedel) ldata.friedel();
    return true;
  }


  template<class T> HKL_data<T>& HKL_data<T>::operator =( const HKL_data<T>& other )
  {
    if ( parent_hkl_info == NULL ) {
      init( other );
    } else {
      if ( parent_hkl_info != other.parent_hkl_info )
        Message::message( Message_fatal( "HKL_data<T>: mismatched parent HKL_info is assignment" ) );
    }
    list = other.list;
    return *this;
  }


  template<class T> HKL_data<T>& HKL_data<T>::operator =( const T& value )
  {
    for ( int i = 0; i < list.size(); i++ ) list[i] = value;
    return *this;
  }


  template<class T> void HKL_data<T>::debug() const
  {
    HKL_data_base::debug();
    std::cout << "Size " << list.size() << "\n";
  }


} // namespace clipper

#endif