00001 #ifndef ARITHMETICPATHCV_H 00002 #define ARITHMETICPATHCV_H 00003 00004 #include "colvarmodule.h" 00005 00006 #include <vector> 00007 #include <cmath> 00008 #include <limits> 00009 #include <string> 00010 00011 namespace ArithmeticPathCV { 00012 00013 using std::vector; 00014 00015 enum path_sz {S, Z}; 00016 00017 template <typename element_type, typename scalar_type, path_sz path_type> 00018 class ArithmeticPathBase { 00019 public: 00020 ArithmeticPathBase() {} 00021 virtual ~ArithmeticPathBase() {} 00022 virtual void initialize(size_t p_num_elements, size_t p_total_frames, double p_lambda, const vector<element_type>& p_element, const vector<double>& p_weights); 00023 virtual void updateDistanceToReferenceFrames() = 0; 00024 virtual void computeValue(); 00025 virtual void computeDerivatives(); 00026 virtual void compute(); 00027 virtual void reComputeLambda(const vector<scalar_type>& rmsd_between_refs); 00028 protected: 00029 scalar_type lambda; 00030 vector<scalar_type> weights; 00031 size_t num_elements; 00032 size_t total_frames; 00033 vector< vector<element_type> > frame_element_distances; 00034 scalar_type s; 00035 scalar_type z; 00036 vector<element_type> dsdx; 00037 vector<element_type> dzdx; 00038 private: 00039 // intermediate variables 00040 vector<scalar_type> s_numerator_frame; 00041 vector<scalar_type> s_denominator_frame; 00042 scalar_type numerator_s; 00043 scalar_type denominator_s; 00044 scalar_type normalization_factor; 00045 }; 00046 00047 template <typename element_type, typename scalar_type, path_sz path_type> 00048 void ArithmeticPathBase<element_type, scalar_type, path_type>::initialize(size_t p_num_elements, size_t p_total_frames, double p_lambda, const vector<element_type>& p_element, const vector<double>& p_weights) { 00049 lambda = p_lambda; 00050 weights = p_weights; 00051 num_elements = p_num_elements; 00052 total_frames = p_total_frames; 00053 frame_element_distances.resize(total_frames, p_element); 00054 for (size_t i_frame = 0; i_frame < frame_element_distances.size(); ++i_frame) { 00055 for (size_t j_elem = 0; j_elem < num_elements; ++j_elem) { 00056 frame_element_distances[i_frame][j_elem].reset(); 00057 } 00058 } 00059 s = scalar_type(0); 00060 z = scalar_type(0); 00061 dsdx = p_element; 00062 dzdx = p_element; 00063 s_numerator_frame.resize(total_frames, scalar_type(0)); 00064 s_denominator_frame.resize(total_frames, scalar_type(0)); 00065 numerator_s = scalar_type(0); 00066 denominator_s = scalar_type(0); 00067 normalization_factor = 1.0 / static_cast<scalar_type>(total_frames - 1); 00068 } 00069 00070 template <typename element_type, typename scalar_type, path_sz path_type> 00071 void ArithmeticPathBase<element_type, scalar_type, path_type>::computeValue() { 00072 updateDistanceToReferenceFrames(); 00073 numerator_s = scalar_type(0); 00074 denominator_s = scalar_type(0); 00075 for (size_t i_frame = 0; i_frame < frame_element_distances.size(); ++i_frame) { 00076 scalar_type exponent_tmp = scalar_type(0); 00077 for (size_t j_elem = 0; j_elem < num_elements; ++j_elem) { 00078 exponent_tmp += weights[j_elem] * frame_element_distances[i_frame][j_elem] * weights[j_elem] * frame_element_distances[i_frame][j_elem]; 00079 } 00080 exponent_tmp = exponent_tmp * -1.0 * lambda; 00081 // prevent underflow if the argument of cvm::exp is less than -708.4 00082 if (exponent_tmp > -708.4) { 00083 exponent_tmp = cvm::exp(exponent_tmp); 00084 } else { 00085 exponent_tmp = 0; 00086 } 00087 numerator_s += static_cast<scalar_type>(i_frame) * exponent_tmp; 00088 denominator_s += exponent_tmp; 00089 s_numerator_frame[i_frame] = static_cast<scalar_type>(i_frame) * exponent_tmp; 00090 s_denominator_frame[i_frame] = exponent_tmp; 00091 } 00092 s = numerator_s / denominator_s * normalization_factor; 00093 z = -1.0 / lambda * cvm::logn(denominator_s); 00094 } 00095 00096 template <typename element_type, typename scalar_type, path_sz path_type> 00097 void ArithmeticPathBase<element_type, scalar_type, path_type>::compute() { 00098 computeValue(); 00099 computeDerivatives(); 00100 } 00101 00102 template <typename element_type, typename scalar_type, path_sz path_type> 00103 void ArithmeticPathBase<element_type, scalar_type, path_type>::computeDerivatives() { 00104 for (size_t j_elem = 0; j_elem < num_elements; ++j_elem) { 00105 element_type dsdxj_numerator_part1(dsdx[j_elem]); 00106 element_type dsdxj_numerator_part2(dsdx[j_elem]); 00107 element_type dzdxj_numerator(dsdx[j_elem]); 00108 dsdxj_numerator_part1.reset(); 00109 dsdxj_numerator_part2.reset(); 00110 dzdxj_numerator.reset(); 00111 for (size_t i_frame = 0; i_frame < frame_element_distances.size(); ++i_frame) { 00112 element_type derivative_tmp = -2.0 * lambda * weights[j_elem] * weights[j_elem] * frame_element_distances[i_frame][j_elem]; 00113 dsdxj_numerator_part1 += s_numerator_frame[i_frame] * derivative_tmp; 00114 dsdxj_numerator_part2 += s_denominator_frame[i_frame] * derivative_tmp; 00115 dzdxj_numerator += s_denominator_frame[i_frame] * derivative_tmp; 00116 } 00117 dsdxj_numerator_part1 *= denominator_s; 00118 dsdxj_numerator_part2 *= numerator_s; 00119 if ((dsdxj_numerator_part1 - dsdxj_numerator_part2).norm() < std::numeric_limits<scalar_type>::min()) { 00120 dsdx[j_elem] = 0; 00121 } else { 00122 dsdx[j_elem] = (dsdxj_numerator_part1 - dsdxj_numerator_part2) / (denominator_s * denominator_s) * normalization_factor; 00123 } 00124 dzdx[j_elem] = -1.0 / lambda * dzdxj_numerator / denominator_s; 00125 } 00126 } 00127 00128 template <typename element_type, typename scalar_type, path_sz path_type> 00129 void ArithmeticPathBase<element_type, scalar_type, path_type>::reComputeLambda(const vector<scalar_type>& rmsd_between_refs) { 00130 scalar_type mean_square_displacements = 0.0; 00131 for (size_t i_frame = 1; i_frame < total_frames; ++i_frame) { 00132 cvm::log(std::string("Distance between frame ") + cvm::to_str(i_frame) + " and " + cvm::to_str(i_frame + 1) + " is " + cvm::to_str(rmsd_between_refs[i_frame - 1]) + std::string("\n")); 00133 mean_square_displacements += rmsd_between_refs[i_frame - 1] * rmsd_between_refs[i_frame - 1]; 00134 } 00135 mean_square_displacements /= scalar_type(total_frames - 1); 00136 lambda = 1.0 / mean_square_displacements; 00137 } 00138 } 00139 00140 #endif // ARITHMETICPATHCV_H