log_scale_spline.h

// -*-c++-*-
#ifndef STORMM_LOG_SCALE_SPLINE_H
#define STORMM_LOG_SCALE_SPLINE_H
 
#include <cmath>
#include <vector>
#include "copyright.h"
#include "Accelerator/gpu_details.h"
#include "Accelerator/hybrid.h"
#include "Constants/behavior.h"
#include "Constants/symbol_values.h"
#include "DataTypes/common_types.h"
#include "DataTypes/stormm_vector_types.h"
#include "Math/matrix_ops.h"
#include "Parsing/parse.h"
#include "Parsing/parsing_enumerators.h"
#include "Potential/pme_util.h"
#include "Reporting/error_format.h"
#include "math_enumerators.h"
 
namespace stormm {
namespace stmath {
 
using card::GpuDetails;
using card::Hybrid;
using card::HybridTargetLevel;
using constants::ExceptionResponse;
using constants::PrecisionModel;
using symbols::amber_ancient_bioq;
using energy::minimum_ewald_coefficient;
using energy::maximum_ewald_coefficient;
using parse::minimalRealFormat;
using parse::NumberFormat;
using parse::realToString;

constexpr int default_logtab_mantissa_bits = 5;
constexpr double default_logtab_max_range = 64.0;
constexpr double default_logtab_min_range = 0.00390625;
  
template <typename T4> struct LogSplineTable {

  LogSplineTable(BasisFunctions basis_in, TableIndexing lookup_in, int detail_bits_in,
                 int index_bound_in, uint sp_detail_mask_in, ullint dp_detail_mask_in,
                 const float arg_offset_in, const T4* table_in);

  LogSplineTable(const LogSplineTable &original) = default;
  LogSplineTable(LogSplineTable &&original) = default;

  const BasisFunctions basis;

  const TableIndexing lookup;
  
  const int detail_bits;

  const int index_bound;

  const uint sp_detail_mask;
  const ullint dp_detail_mask;

  const float arg_offset;
  
  const T4* table;
};
  
template <typename T4> class LogScaleSpline {
public:

  LogScaleSpline(LogSplineForm target_form_in, double ewald_coefficient_in,
                 double coulomb_constant_in = amber_ancient_bioq,
                 int mantissa_bits_in = default_logtab_mantissa_bits,
                 double max_range_in = default_logtab_max_range,
                 double min_range_in = default_logtab_min_range,
                 TableIndexing indexing_method_in = TableIndexing::SQUARED_ARG,
                 BasisFunctions basis_set_in = BasisFunctions::POLYNOMIAL,
                 int ulp_optimization_depth_in = 2,
                 float indexing_offset_in = 0.0,
                 ExceptionResponse policy = ExceptionResponse::DIE);
 
  LogScaleSpline(LogSplineForm target_form_in,
                 const std::vector<std::vector<double2>> &custom_form_in,
                 int mantissa_bits_in = default_logtab_mantissa_bits,
                 double max_range_in = default_logtab_max_range,
                 double min_range_in = default_logtab_min_range,
                 TableIndexing indexing_method_in = TableIndexing::SQUARED_ARG,
                 BasisFunctions basis_set_in = BasisFunctions::POLYNOMIAL,
                 int ulp_optimization_depth_in = 2,
                 float indexing_offset_in = 0.0,
                 ExceptionResponse policy = ExceptionResponse::WARN);

  LogScaleSpline(const LogScaleSpline &original) = default;
  LogScaleSpline(LogScaleSpline &&original) = default;
  LogScaleSpline& operator=(const LogScaleSpline &original) = default;
  LogScaleSpline& operator=(LogScaleSpline &&original) = default;

  LogSplineForm getForm() const;

  TableIndexing getIndexingMethod() const;

  BasisFunctions getBasisSet() const;
  
  int getBitStride() const;

  int getSplineDensity() const;

  double getEwaldCoefficient() const;

  double getCoulombConstant() const;
  
  double getMaximumRange() const;

  double getMinimumRange() const;

  float getIndexingOffset() const;

  int getOptimizationDepth() const;
  
  int getSplineIndex(double r) const;

  int getSplineIndexByRealArg(double r) const;

  double3 getErrorEstimate() const;
  double3 getErrorEstimate(double r) const;

  double interpolateByValue(double r) const;

  double interpolateBySquaredValue(double r2) const;

  double evaluate(double r) const;

  double evaluateByRealArg(double r) const;  
  
  const LogSplineTable<T4> data(HybridTargetLevel tier = HybridTargetLevel::HOST) const;

  const LogSplineTable<void>
  templateFreeData(HybridTargetLevel tier = HybridTargetLevel::HOST) const;

  const LogScaleSpline<T4>* getSelfPointer() const;
  
#ifdef STORMM_USE_HPC
  void upload();
  
  void download();
#endif
  void setMinimumSignificantRange(double min_range_in);
 
private:
  LogSplineForm target_form;                
  TableIndexing indexing_method;            
  BasisFunctions basis_set;                 
  int mantissa_bits;                        
  ullint dp_detail_bitmask;                 
  uint sp_detail_bitmask;                   
  int segments_per_stride;                  
  double ewald_coefficient;                 
  double coulomb_constant;                  
  double maximum_range;                     
  double minimum_absolute_range;            
  double minimum_significant_range;         
  float indexing_offset;                    
  int ulp_optimization_depth;               
  PrecisionModel precision;                 
  Hybrid<T4> table;                         
  double mean_overall_error;                
  double stdev_overall_error;               
  double max_overall_error;                 
  std::vector<double> mean_segment_error;   
  std::vector<double> stdev_segment_error;  
  std::vector<double> max_segment_error;    

  void setPrecisionModel();
  
  void setTargetForm(LogSplineForm target_form_in,
                     const std::vector<std::vector<double2>> &custom_form_in = {});

  void setEwaldCoefficient(double ewald_coefficient_in,
                           ExceptionResponse policy = ExceptionResponse::DIE);

  void setMantissaBits(int mantissa_bits_in,
                       ExceptionResponse policy = ExceptionResponse::WARN);

  void setMinimumAbsoluteRange(double min_range_in);

  void setULPDepth(int ulp_optimization_depth_in);
  
  float checkIndexingOffset(float indexing_offset_in);
  
  void allocate(double max_range_in, const std::vector<std::vector<double2>> &custom_form_in = {},
                const ExceptionResponse policy = ExceptionResponse::WARN);

  void evaluatePrescribedTable();

  void evaluateCustomTable();

  void evaluateTargetInSegment(std::vector<double> *bwrk, std::vector<double> *bref, double r2_i,
                               double width);

  void evaluateSplineInSegment(std::vector<double> *estm, const T4 spl_tmp, double r2_i,
                               double width);

  void evaluateOverallError();
};

ullint doublePrecisionSplineDetailMask(int mantissa_bits);

uint singlePrecisionSplineDetailMask(int mantissa_bits);

double doublePrecisionSplineArgument(double arg, int mantissa_bits, ullint detail_mask);

float singlePrecisionSplineArgument(float arg, int mantissa_bits, uint detail_mask);

template <typename T4> LogSplineTable<T4> restoreType(const LogSplineTable<void> *rasa);
 
template <typename T4> LogSplineTable<T4> restoreType(const LogSplineTable<void> &rasa);
  
} // namespace stmath
} // namespace stormm
 
#include "log_scale_spline.tpp"
 
#endif