synthesis_permutor.h

// -*-c++-*-
#ifndef STORMM_SYNTHESIS_PERMUTOR_H
#define STORMM_SYNTHESIS_PERMUTOR_H
 
#include <vector>
#include "copyright.h"
#include "Accelerator/hybrid.h"
#include "Chemistry/chemical_features.h"
#include "Chemistry/chemistry_enumerators.h"
#include "Constants/fixed_precision.h"
#include "Math/tickcounter.h"
#include "Namelists/nml_conformer.h"
#include "Random/random.h"
#include "Structure/clash_detection.h"
#include "Structure/isomerization.h"
#include "Structure/local_arrangement.h"
#include "Structure/structure_enumerators.h"
#include "Topology/atomgraph.h"
#include "atomgraph_synthesis.h"
#include "phasespace_synthesis.h"
#include "synthesis_enumerators.h"
 
namespace stormm {
namespace synthesis {
 
using card::Hybrid;
using card::HybridTargetLevel;
using chemistry::ChemicalFeatures;
using chemistry::ChiralInversionProtocol;
using chemistry::ChiralOrientation;
using chemistry::getChiralOrientation;
using chemistry::ConformationEdit;
using chemistry::CoupledEdit;
using chemistry::IsomerPlan;
using namelist::ConformerControls;
using namelist::default_conf_clash_pairs;
using namelist::default_conf_max_seeding_attempts;
using numerics::default_globalpos_scale_bits;
using numerics::default_localpos_scale_bits;
using numerics::default_velocity_scale_bits;
using numerics::default_force_scale_bits;
using random::Xoshiro256ppGenerator;
using stmath::TickCounter;
using structure::ClashReport;
using structure::detectClash;
using structure::dihedralAngle;
using structure::flipChiralCenter;
using structure::rotateAboutBond;
using structure::SamplingIntensity;
using topology::AtomGraph;

template <typename T> struct SyPermutorKit {
public:

  explicit SyPermutorKit(int system_count_in, int perm_map_count_in,
                         const int* perm_map_idx_in, const int2* perm_elements_in,
                         const int* perm_element_bounds_in, const int* system_settings_in,
                         const int* system_settings_limits_in, const int* rot_grp_atoms_in,
                         const int* rot_grp_bounds_in, const int* prm_rot_grp_bounds_in,
                         const int* ctx_grp_atoms_in, const int* ctx_grp_bounds_in,
                         const int* prm_ctx_grp_bounds_in, const int* inv_grp_atoms_in,
                         const int* inv_grp_bounds_in, const int* prm_inv_grp_bounds_in,
                         const int* chiral_atoms_in, const int* chiral_protocols_in,
                         const int4* rot_bond_markers_in, const int4* ctx_bond_markers_in,
                         const int4* chiral_markers_in, const T* rot_bond_settings_in,
                         const T* ctx_bond_settings_in, const int* rot_bond_settings_bounds_in,
                         const int* ctx_bond_settings_bounds_in,
                         const int* chiral_settings_in, const int* chiral_settings_bounds_in);
 
  // The member variables store size constants for arrays (the numbers of systems in the synthesis
  // and the number of unique prmologies (translating to the number of unique permutor maps)
  // involved in manipulating each structure of the synthesis.
  int system_count;                     
  int perm_map_count;                   
  const int* perm_map_idx;              
  const int2* perm_elements;            
  const int* perm_element_bounds;       
  const int* system_settings;           
  const int* system_settings_limits;    
  const int* rot_grp_atoms;             
  const int* rot_grp_bounds;            
  const int* prm_rot_grp_bounds;        
  const int* ctx_grp_atoms;             
  const int* ctx_grp_bounds;            
  const int* prm_ctx_grp_bounds;        
  const int* inv_grp_atoms;             
  const int* inv_grp_bounds;            
  const int* prm_inv_grp_bounds;        
  const int* chiral_atoms;              
  const int* chiral_protocols;          
  const int4* rot_bond_markers;         
  const int4* ctx_bond_markers;         
  const int4* chiral_markers;           
  const T* rot_bond_settings;           
  const T* ctx_bond_settings;           
  const int* chiral_settings;           
  const int* rot_bond_settings_bounds;  
  const int* ctx_bond_settings_bounds;  
  const int* chiral_settings_bounds;    
};
  
class SynthesisPermutor {
public:

  SynthesisPermutor();
 
  SynthesisPermutor(const AtomGraphSynthesis *poly_ag, StopWatch *timer = nullptr);
 
  SynthesisPermutor(const AtomGraphSynthesis &poly_ag, StopWatch *timer = nullptr);
 
  SynthesisPermutor(const AtomGraphSynthesis *poly_ag,
                    const ConformerControls &confcon, StopWatch *timer = nullptr);
 
  SynthesisPermutor(const AtomGraphSynthesis &poly_ag,
                    const ConformerControls &confcon, StopWatch *timer = nullptr);
  
  SynthesisPermutor(const PhaseSpaceSynthesis *poly_ps, StopWatch *timer = nullptr);
 
  SynthesisPermutor(const PhaseSpaceSynthesis &poly_ps, StopWatch *timer = nullptr);
 
  SynthesisPermutor(const PhaseSpaceSynthesis *poly_ps,
                    const ConformerControls &confcon, StopWatch *timer = nullptr);
 
  SynthesisPermutor(const PhaseSpaceSynthesis &poly_ps,
                    const ConformerControls &confcon, StopWatch *timer = nullptr);
  
  SynthesisPermutor(const std::vector<ChemicalFeatures*> &chemfe_in, bool retain_pointers = true,
                    StopWatch *timer = nullptr);
 
  SynthesisPermutor(const std::vector<ChemicalFeatures> &chemfe_in, bool retain_pointers = true,
                    StopWatch *timer = nullptr);
 
  SynthesisPermutor(const std::vector<ChemicalFeatures*> &chemfe_in,
                    const PhaseSpaceSynthesis *poly_ps, const ConformerControls &confcon,
                    const bool retain_pointers = true, StopWatch *timer = nullptr);
 
  SynthesisPermutor(const std::vector<ChemicalFeatures*> &chemfe_in,
                    const PhaseSpaceSynthesis &poly_ps, const ConformerControls &confcon,
                    const bool retain_pointers = true, StopWatch *timer = nullptr);
 
  SynthesisPermutor(const std::vector<ChemicalFeatures> &chemfe_in,
                    const PhaseSpaceSynthesis &poly_ps, const ConformerControls &confcon,
                    const bool retain_pointers = true, StopWatch *timer = nullptr);

  SynthesisPermutor(const SynthesisPermutor &original);
  SynthesisPermutor(SynthesisPermutor &&original);
  SynthesisPermutor& operator=(const SynthesisPermutor &other);
  SynthesisPermutor& operator=(SynthesisPermutor &&other);
  
  int getSystemCount() const;

  int getPermutorSetCount() const;

  int getPermutorMapIndex(const AtomGraph *query_ag) const;
  int getPermutorMapIndex(const AtomGraph &query_ag) const;
  int getPermutorMapIndex(const PhaseSpaceSynthesis *poly_ps, int system_index) const;
  int getPermutorMapIndex(const PhaseSpaceSynthesis &poly_ps, int system_index) const;
  
  int getRotatableBondSampleCount() const;
  
  int getCisTransBondSampleCount() const;

  int getRotatableBondCount(const int permutor_map_index) const;
  int getRotatableBondCount() const;

  int getCisTransBondCount(const int permutor_map_index) const;
  int getCisTransBondCount() const;
  
  int getChiralCenterCount(const int permutor_map_index) const;
  int getChiralCenterCount() const;
  
  const std::vector<int>& getElementSampleCounts(int system_index) const;
  const std::vector<int>& getElementSampleCounts(const AtomGraph *query_ag) const;
  const std::vector<int>& getElementSampleCounts(const AtomGraph &query_ag) const;

  int getVariableElementCount(int system_index) const;

  int getElementSampleCount(ConformationEdit element_kind, int index) const;
  int getElementSampleCount(CoupledEdit ce) const;

  int getSystemRotatableBondCount(int system_index) const;

  int getSystemCisTransBondCount(int system_index) const;

  int getSystemChiralCenterCount(int system_index) const;

  const IsomerPlan& getRotatableGroup(int system_index, int group_index) const;

  const IsomerPlan& getCisTransGroup(int system_index, int group_index) const;

  const IsomerPlan& getInvertibleGroup(int system_index, int group_index) const;

  const TickCounter<double>& getStateTracker(int system_index) const;
  const TickCounter<double>& getStateTracker(const AtomGraph *ag) const;
  const TickCounter<double>& getStateTracker(const AtomGraph &ag) const;

  llint getReplicaCount(int query_index, SamplingIntensity effort) const;
  llint getReplicaCount(const AtomGraph *query_ag, SamplingIntensity effort) const;
  llint getReplicaCount(const AtomGraph &query_ag, SamplingIntensity effort) const;
  std::vector<llint> getReplicaCount(SamplingIntensity effort) const;

  const ChemicalFeatures* getChemicalFeaturesPointer(int map_index) const;
  const std::vector<ChemicalFeatures*> getChemicalFeaturesPointer() const;
  
  const PhaseSpaceSynthesis* getSynthesisPointer() const;
  
  const SyPermutorKit<double> dpData(HybridTargetLevel tier = HybridTargetLevel::HOST) const;

  const SyPermutorKit<float> spData(HybridTargetLevel tier = HybridTargetLevel::HOST) const;
 
#ifdef STORMM_USE_HPC
  void upload();

  void download();
#endif

  void defineRotatableBondSettings(int map_index,
                                   const std::vector<std::vector<double>> &settings);

  void defineCisTransBondSettings(int system_index,
                                  const std::vector<std::vector<double>> &settings);

  void defineChiralCenterSettings(int system_index,
                                  const std::vector<std::vector<int>> &settings);

  void applySynthesis(const PhaseSpaceSynthesis *poly_ps_in,
                      VariableTorsionAdjustment adj = VariableTorsionAdjustment::DO_NOT_CHANGE);
 
  void applySynthesis(const PhaseSpaceSynthesis &poly_ps_in,
                      VariableTorsionAdjustment adj = VariableTorsionAdjustment::DO_NOT_CHANGE);
  
  void setVariableRanges(const std::vector<double> &general_rot_settings_in,
                         const std::vector<double> &general_ctx_settings_in,
                         bool sample_cis_trans_states = false, bool sample_chiral_states = false);

  void impartControlData(const ConformerControls &confcon);

  void permuteSystem(PhaseSpaceSynthesis *psynth, int system_index, int map_index,
                     const CoupledEdit ce, int setting_index,
                     PrecisionModel prec = PrecisionModel::SINGLE) const;
 
  void permuteSystem(int system_index, const CoupledEdit ce, int setting_index,
                     PrecisionModel prec = PrecisionModel::SINGLE);

  void randomizeSystem(PhaseSpaceSynthesis *psynth, int system_index, int map_index,
                       Xoshiro256ppGenerator *xrs, PrecisionModel prec = PrecisionModel::SINGLE,
                       const CoupledEdit mask_a = { ConformationEdit::BOND_ROTATION, -1 },
                       const CoupledEdit mask_b = { ConformationEdit::BOND_ROTATION, -1 },
                       const CoupledEdit mask_c = { ConformationEdit::BOND_ROTATION, -1 }) const;
  
  void randomizeSystem(int system_index, Xoshiro256ppGenerator *xrs,
                       PrecisionModel prec = PrecisionModel::SINGLE,
                       const CoupledEdit mask_a = { ConformationEdit::BOND_ROTATION, -1 },
                       const CoupledEdit mask_b = { ConformationEdit::BOND_ROTATION, -1 },
                       const CoupledEdit mask_c = { ConformationEdit::BOND_ROTATION, -1 });

  PhaseSpaceSynthesis buildSynthesis(SamplingIntensity effort, Xoshiro256ppGenerator *xrs,
                                     std::vector<int> *correspondence, int system_limit,
                                     int globalpos_scale_bits = default_globalpos_scale_bits,
                                     int localpos_scale_bits = default_localpos_scale_bits,
                                     int velocity_scale_bits = default_velocity_scale_bits,
                                     int force_scale_bits = default_force_scale_bits,
                                     PrecisionModel prec = PrecisionModel::SINGLE,
                                     ClashReport *clrep = nullptr,
                                     int iteration_limit = default_conf_max_seeding_attempts,
                                     int max_clashes = default_conf_clash_pairs) const;
 
  PhaseSpaceSynthesis buildSynthesis(SamplingIntensity effort, Xoshiro256ppGenerator *xrs,
                                     std::vector<int> *correspondence, PrecisionModel prec,
                                     int system_limit, ClashReport *clrep = nullptr,
                                     int iteration_limit = default_conf_max_seeding_attempts,
                                     int max_clashes = default_conf_clash_pairs) const;
 
  PhaseSpaceSynthesis buildSynthesis(SamplingIntensity effort, Xoshiro256ppGenerator *xrs,
                                     int system_limit,
                                     int globalpos_scale_bits = default_globalpos_scale_bits,
                                     int localpos_scale_bits = default_localpos_scale_bits,
                                     int velocity_scale_bits = default_velocity_scale_bits,
                                     int force_scale_bits = default_force_scale_bits,
                                     PrecisionModel prec = PrecisionModel::SINGLE,
                                     ClashReport *clrep = nullptr,
                                     int iteration_limit = default_conf_max_seeding_attempts,
                                     int max_clashes = default_conf_clash_pairs) const;
 
  PhaseSpaceSynthesis buildSynthesis(SamplingIntensity effort, Xoshiro256ppGenerator *xrs,
                                     PrecisionModel prec, int system_limit,
                                     ClashReport *clrep = nullptr,
                                     int iteration_limit = default_conf_max_seeding_attempts,
                                     int max_clashes = default_conf_clash_pairs) const;
 
private:
  int system_count;                              
  int permutor_map_count;                        
  int rotatable_bond_samples;                    
  int cis_trans_bond_samples;                    
  bool features_objects_accessible;              
  double rotation_setting_snap;                  
  double cis_trans_setting_snap;                 
  int total_rotatable_groups;                    
  int total_cis_trans_groups;                    
  int total_invertible_groups;                   
  Hybrid<int> permutor_map_indices;              
  Hybrid<int2> permutor_elements;                
  Hybrid<int> permutor_element_bounds;           
  Hybrid<int> system_settings;                   
  Hybrid<int> system_settings_limits;            
  Hybrid<int> synthesis_data;                    
  Hybrid<int> rotatable_group_atoms;             
  Hybrid<int> rotatable_group_bounds;            
  Hybrid<int> permutor_rotatable_group_bounds;   
  Hybrid<int> cis_trans_group_atoms;             
  Hybrid<int> cis_trans_group_bounds;            
  Hybrid<int> permutor_cis_trans_group_bounds;   
  Hybrid<int> invertible_group_atoms;            
  Hybrid<int> invertible_group_bounds;           
  Hybrid<int> permutor_invertible_group_bounds;  
  Hybrid<int> invertible_atom_centers;           
  Hybrid<int> invertible_group_protocols;        
  Hybrid<int4> rotatable_bond_markers;           
  Hybrid<int4> cis_trans_bond_markers;           
  Hybrid<int4> chiral_markers;                   
  Hybrid<double> rotatable_bond_settings;        
  Hybrid<double> cis_trans_bond_settings;        
  Hybrid<float> sp_rotatable_bond_settings;      
  Hybrid<float> sp_cis_trans_bond_settings;      
  Hybrid<int> rotatable_bond_settings_bounds;    
  Hybrid<int> cis_trans_bond_settings_bounds;    
  Hybrid<int> chiral_settings_bounds;            
  Hybrid<int> group_data;                        
  Hybrid<int4> marker_data;                      
  Hybrid<int> chiral_settings;                   
 
  // The following arrays store couple variables needed to understand "light" and "heavy" sampling.
  // The case of "minimal" sampling is trivial--sample all states of each variable individually
  // while other variables get set to random values.  In "light" sampling the goal is to find all
  // combinations of two coupled variables, but some variables will exist only in isolation.  In
  // "heavy" sampling combinations of three variables are applied, but some variables may have no
  // couplings or only one other variable to connect with.  For these cases where pairs or triplets
  // cannot be fulfilled, the tuples representing them will carry a value of -1.  Exhaustive
  // sampling will place all variables in one massive tuple and do all permutations.
  Hybrid<int2> light_sampling_foci;        
  Hybrid<int2> heavy_sampling_foci;        
  Hybrid<int> light_sampling_foci_bounds;  
  Hybrid<int> heavy_sampling_foci_bounds;  
  Hybrid<int2> variables_data;             
  
  // The following CPU-based data arrays help to stage the Hybrid data.
  std::vector<double> general_rotation_settings;  
  std::vector<double> general_cis_trans_settings; 
  std::vector<IsomerPlan> rotatable_groups;       
  std::vector<IsomerPlan> cis_trans_groups;       
  std::vector<IsomerPlan> invertible_groups;      
 
  // These CPU-based arrays put dimensions on the permutations required to explore each system.
  std::vector<TickCounter<double>> state_trackers;  
  std::vector<int> minimal_sampling_replicas;       
  std::vector<int> light_sampling_replicas;         
  std::vector<int> heavy_sampling_replicas;         
  std::vector<llint> exhaustive_sampling_replicas;  

  std::vector<AtomGraph*> topologies;
  
  std::vector<ChemicalFeatures*> features;       

  PhaseSpaceSynthesis *poly_ps_ptr;
  
  void rebasePointers();
  
  std::vector<ChemicalFeatures> temporaryFeatures(StopWatch *timer = nullptr) const;

  const std::vector<ChemicalFeatures*>
  temporaryFeaturesPointers(const std::vector<ChemicalFeatures> &feat_in) const;

  void fillPermutorMaps(const std::vector<ChemicalFeatures*> &chemfe_in, StopWatch *timer);
  
  void validateSystemIndex(int system_index) const;

  void validateMapIndex(int map_index) const;

  template <typename T> void alterMapSettings(int map_index, Hybrid<T> *variable_settings,
                                              const Hybrid<int> &variable_settings_bounds,
                                              const Hybrid<int> &permutor_map_bounds,
                                              const std::vector<std::vector<T>> &settings);

  template <typename T>
  void permuteSystem(CoordinateSeries<T> *cs, int map_index, CoupledEdit ce,
                     int setting_index) const;

  template <typename T>
  void randomizeSystem(CoordinateSeries<T> *cs, int map_index, Xoshiro256ppGenerator *xrs,
                       const CoupledEdit mask_a = { ConformationEdit::BOND_ROTATION, -1 },
                       const CoupledEdit mask_b = { ConformationEdit::BOND_ROTATION, -1 },
                       const CoupledEdit mask_c = { ConformationEdit::BOND_ROTATION, -1 }) const;

  template <typename T>
  bool resolveClashes(CoordinateSeries<T> *cs, int map_index, Xoshiro256ppGenerator *xrs,
                      const StaticExclusionMask &excl, ClashReport *summary, int iteration_limit,
                      int max_clashes,
                      const CoupledEdit mask_a = { ConformationEdit::BOND_ROTATION, -1 },
                      const CoupledEdit mask_b = { ConformationEdit::BOND_ROTATION, -1 },
                      const CoupledEdit mask_c = { ConformationEdit::BOND_ROTATION, -1 }) const;
};
 
} // namespace synthesis
} // namespace stormm
 
#include "synthesis_permutor.tpp"
 
#endif