Main Page   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members   Related Pages  

colvar_arithmeticpath.h

Go to the documentation of this file.
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

Generated on Wed Oct 9 02:42:05 2024 for VMD (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002