Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | File List | Class Members | File Members

colvargrid.h

Go to the documentation of this file.
00001 // -*- c++ -*-
00002 
00003 #ifndef COLVARGRID_H
00004 #define COLVARGRID_H
00005 
00006 #include <iostream>
00007 #include <iomanip>
00008 #include <cmath>
00009 
00010 #include "colvar.h"
00011 #include "colvarmodule.h"
00012 #include "colvarvalue.h"
00013 #include "colvarparse.h"
00014 
00019 template <class T> class colvar_grid : public colvarparse {
00020 
00021 protected:
00022 
00024   size_t nd;
00025 
00027   std::vector<int> nx;
00028 
00030   std::vector<int> nxc;
00031 
00034   size_t mult;
00035 
00037   size_t nt;
00038 
00040   std::vector<T> data;
00041 
00043   std::vector<size_t> new_data;
00044 
00046   std::vector<colvar *> cv;
00047 
00049   inline size_t address (std::vector<int> const &ix) const
00050   {
00051     size_t addr = 0;
00052     for (size_t i = 0; i < nd; i++) {
00053       addr += ix[i]*nxc[i];
00054       if (cvm::debug()) {
00055         if (ix[i] >= nx[i])
00056           cvm::fatal_error ("Error: exceeding bounds in colvar_grid.\n");
00057       }
00058     }
00059     return addr;
00060   }
00061 
00062 public:
00063 
00065   std::vector<colvarvalue> lower_boundaries;
00066 
00068   std::vector<colvarvalue> upper_boundaries;
00069 
00071   std::vector<bool>        periodic;
00072 
00074   std::vector<cvm::real>   widths;
00075 
00077   bool has_parent_data;
00078 
00080   bool has_data;
00081 
00083   inline size_t number_of_colvars() const
00084   {
00085     return nd;
00086   }
00087 
00090   inline size_t number_of_points (int const icv = -1) const
00091   {
00092     if (icv < 0) {
00093       return nt;
00094     } else {
00095       return nx[icv];
00096     }
00097   }
00098 
00100   inline std::vector<int> const & sizes() const
00101   {
00102     return nx;
00103   }
00104 
00106   inline void set_sizes (std::vector<int> const &new_sizes)
00107   {
00108     nx = new_sizes;
00109   }
00110 
00112   inline size_t multiplicity() const
00113   {
00114     return mult;
00115   }
00116 
00118   void create (std::vector<int> const &nx_i,
00119                T const &t = T(),
00120                size_t const &mult_i = 1)
00121   {
00122     mult = mult_i;
00123     nd = nx_i.size();
00124     nxc.resize (nd);
00125     nx = nx_i;
00126 
00127     nt = mult;
00128     for (int i = nd-1; i >= 0; i--) {
00129       if (nx_i[i] <= 0)
00130         cvm::fatal_error ("Error: providing an invalid number of points, "+
00131                           cvm::to_str (nx_i[i])+".\n");
00132       nxc[i] = nt;
00133       nt *= nx[i];
00134     }
00135 
00136     data.reserve (nt);
00137     data.assign (nt, t);
00138   }
00139 
00141   void create()
00142   {
00143     create (this->nx, T(), this->mult);
00144   }
00145 
00147   void reset (T const &t = T())
00148   {
00149     data.assign (nt, t);
00150   }
00151 
00152 
00154   colvar_grid() : has_data (false)
00155   {
00156     save_delimiters = false;
00157     nd = nt = 0;
00158   }
00159 
00161   virtual ~colvar_grid()
00162   {}
00163 
00167   colvar_grid (colvar_grid<T> const &g) : has_data (false),
00168                                           nd (g.nd),
00169                                           mult (g.mult),
00170                                           cv (g.cv),
00171                                           lower_boundaries (g.lower_boundaries),
00172                                           upper_boundaries (g.upper_boundaries),
00173                                           periodic (g.periodic),
00174                                           widths (g.widths),
00175                                           data()
00176   {
00177     save_delimiters = false;
00178   }
00179 
00184   colvar_grid (std::vector<int> const &nx_i,
00185                T const &t = T(),
00186                size_t const &mult_i = 1) : has_data (false)
00187   {
00188     save_delimiters = false;
00189     this->create (nx_i, t, mult_i);
00190   }
00191 
00193   colvar_grid (std::vector<colvar *> const &colvars,
00194                T const &t = T(),
00195                size_t const &mult_i = 1,
00196                bool margin = false)
00197     : cv (colvars), has_data (false)
00198   {
00199     save_delimiters = false;
00200 
00201     std::vector<int> nx_i;
00202 
00203     if (cvm::debug())
00204       cvm::log ("Allocating a grid for "+cvm::to_str (colvars.size())+
00205                 " collective variables.\n");
00206 
00207     for (size_t i =  0; i < cv.size(); i++) {
00208 
00209       if (cv[i]->type() != colvarvalue::type_scalar) {
00210         cvm::fatal_error ("Colvar grids can only be automatically "
00211                           "constructed for scalar variables.  "
00212                           "ABF and histogram can not be used; metadynamics "
00213                           "can be used with useGrids disabled.\n");
00214       }
00215 
00216       if (cv[i]->width <= 0.0) {
00217         cvm::fatal_error ("Tried to initialize a grid on a "
00218                           "variable with negative or zero width.\n");
00219       }
00220 
00221       if (!cv[i]->tasks[colvar::task_lower_boundary] || !cv[i]->tasks[colvar::task_upper_boundary]) {
00222         cvm::fatal_error ("Tried to initialize a grid on a "
00223                           "variable with undefined boundaries.\n");
00224       }
00225 
00226       widths.push_back (cv[i]->width);
00227       periodic.push_back (cv[i]->periodic_boundaries());
00228 
00229       if (margin) {
00230         if (periodic[i]) {
00231           // Shift the grid by half the bin width (values at edges instead of center of bins)
00232           lower_boundaries.push_back (cv[i]->lower_boundary.real_value - 0.5 * widths[i]);
00233           upper_boundaries.push_back (cv[i]->upper_boundary.real_value - 0.5 * widths[i]);
00234         } else {
00235           // Make this grid larger by one bin width
00236           lower_boundaries.push_back (cv[i]->lower_boundary.real_value - 0.5 * widths[i]);
00237           upper_boundaries.push_back (cv[i]->upper_boundary.real_value + 0.5 * widths[i]);
00238         }
00239       } else {
00240         lower_boundaries.push_back (cv[i]->lower_boundary);
00241         upper_boundaries.push_back (cv[i]->upper_boundary);
00242       }
00243 
00244 
00245       { 
00246         cvm::real nbins = ( upper_boundaries[i].real_value -
00247                             lower_boundaries[i].real_value ) / widths[i];
00248         int nbins_round = (int)(nbins+0.5);
00249         
00250         if (std::fabs (nbins - cvm::real (nbins_round)) > 1.0E-10) {
00251           cvm::log ("Warning: grid interval ("+
00252                     cvm::to_str (lower_boundaries[i], cvm::cv_width, cvm::cv_prec)+" - "+
00253                     cvm::to_str (upper_boundaries[i], cvm::cv_width, cvm::cv_prec)+
00254                     ") is not commensurate to its bin width ("+
00255                     cvm::to_str (widths[i], cvm::cv_width, cvm::cv_prec)+").\n");
00256           upper_boundaries[i].real_value = lower_boundaries[i].real_value +
00257             (nbins_round * widths[i]);
00258         }
00259 
00260         if (cvm::debug())
00261           cvm::log ("Number of points is "+cvm::to_str ((int) nbins_round)+
00262                     " for the colvar no. "+cvm::to_str (i+1)+".\n");
00263         
00264         nx_i.push_back (nbins_round);
00265       }
00266 
00267     }
00268 
00269     create (nx_i, t, mult_i);
00270   }
00271 
00272 
00275   inline void wrap (std::vector<int> & ix) const
00276   {
00277     for (size_t i = 0; i < nd; i++) {
00278       if (periodic[i]) {
00279         ix[i] = (ix[i] + nx[i]) % nx[i]; //to ensure non-negative result
00280       } else {
00281         if (ix[i] < 0 || ix[i] >= nx[i])
00282           cvm::fatal_error ("Trying to wrap illegal index vector (non-PBC): "
00283                             + cvm::to_str (ix));
00284       }
00285     }
00286   }
00287 
00289   inline int current_bin_scalar(int const i) const
00290   {
00291     return value_to_bin_scalar (cv[i]->value(), i);
00292   }
00293 
00296   inline int value_to_bin_scalar (colvarvalue const &value, const int i) const
00297   {
00298     return (int) std::floor ( (value.real_value - lower_boundaries[i].real_value) / widths[i] );
00299   }
00300 
00302   inline int value_to_bin_scalar (colvarvalue const &value,
00303                                   colvarvalue const &new_offset,
00304                                   cvm::real   const &new_width) const
00305   {
00306     return (int) std::floor ( (value.real_value - new_offset.real_value) / new_width );
00307   }
00308 
00311   inline colvarvalue bin_to_value_scalar (int const &i_bin, int const i) const
00312   {
00313     return lower_boundaries[i].real_value + widths[i] * (0.5 + i_bin);
00314   }
00315   
00317   inline colvarvalue bin_to_value_scalar (int const &i_bin,
00318                                           colvarvalue const &new_offset,
00319                                           cvm::real const &new_width) const
00320   {
00321     return new_offset.real_value + new_width * (0.5 + i_bin);
00322   }
00323 
00325   inline void set_value (std::vector<int> const &ix,
00326                          T const &t,
00327                          size_t const &imult = 0)
00328   { 
00329     data[this->address (ix)+imult] = t;
00330     has_data = true;
00331   }
00332 
00333 
00336   inline T const & value (std::vector<int> const &ix,
00337                           size_t const &imult = 0) const
00338   {
00339     return data[this->address (ix) + imult];
00340   }
00341 
00342 
00344   inline void add_constant (T const &t)
00345   {
00346     for (size_t i = 0; i < nt; i++) 
00347       data[i] += t;
00348     has_data = true;
00349   }
00350 
00352   inline void multiply_constant (cvm::real const &a)
00353   {
00354     for (size_t i = 0; i < nt; i++) 
00355       data[i] *= a;
00356   }
00357 
00358  
00361   inline std::vector<int> const get_colvars_index (std::vector<colvarvalue> const &values) const
00362   {
00363     std::vector<int> index = new_index();
00364     for (size_t i = 0; i < nd; i++) {
00365       index[i] = value_to_bin_scalar (values[i], i);
00366     }
00367     return index;
00368   }
00369 
00372   inline std::vector<int> const get_colvars_index() const
00373   {
00374     std::vector<int> index = new_index();
00375     for (size_t i = 0; i < nd; i++) {
00376       index[i] = current_bin_scalar (i);
00377     }
00378     return index;
00379   }
00380 
00384   inline cvm::real bin_distance_from_boundaries (std::vector<colvarvalue> const &values)
00385   {
00386     cvm::real minimum = 1.0E+16;
00387     for (size_t i = 0; i < nd; i++) {
00388 
00389       if (periodic[i]) continue;
00390 
00391       cvm::real dl = std::sqrt (cv[i]->dist2 (values[i], lower_boundaries[i])) / widths[i];
00392       cvm::real du = std::sqrt (cv[i]->dist2 (values[i], upper_boundaries[i])) / widths[i];
00393 
00394       if (values[i].real_value < lower_boundaries[i])
00395         dl *= -1.0;
00396       if (values[i].real_value > upper_boundaries[i])
00397         du *= -1.0;
00398 
00399       if (dl < minimum) 
00400         minimum = dl;
00401       if (du < minimum) 
00402         minimum = du;
00403     }
00404 
00405     return minimum;
00406   }
00407 
00408 
00413   void map_grid (colvar_grid<T> const &other_grid)
00414   {
00415     if (other_grid.multiplicity() != this->multiplicity())
00416       cvm::fatal_error ("Error: trying to merge two grids with values of "
00417                         "different multiplicity.\n");
00418 
00419     std::vector<colvarvalue> const &gb  = this->lower_boundaries;
00420     std::vector<cvm::real> const &gw    = this->widths;
00421     std::vector<colvarvalue> const &ogb = other_grid.lower_boundaries;
00422     std::vector<cvm::real> const &ogw   = other_grid.widths;
00423 
00424     std::vector<int> ix = this->new_index();
00425     std::vector<int> oix = other_grid.new_index();
00426 
00427     if (cvm::debug())
00428       cvm::log ("Remapping grid...\n");
00429     for ( ; this->index_ok (ix); this->incr (ix)) {
00430 
00431       for (size_t i = 0; i < nd; i++) {
00432         oix[i] =
00433           value_to_bin_scalar (bin_to_value_scalar (ix[i], gb[i], gw[i]),
00434                                ogb[i],
00435                                ogw[i]);
00436       }
00437 
00438       if (! other_grid.index_ok (oix)) {
00439         continue;
00440       }
00441 
00442       for (size_t im = 0; im < mult; im++) {
00443         this->set_value (ix, other_grid.value (oix, im), im);
00444       }
00445     }
00446 
00447     has_data = true;
00448     if (cvm::debug())
00449       cvm::log ("Remapping done.\n");
00450   }
00451 
00454   void add_grid (colvar_grid<T> const &other_grid,
00455                  cvm::real scale_factor = 1.0)
00456   {
00457     if (other_grid.multiplicity() != this->multiplicity())
00458       cvm::fatal_error ("Error: trying to sum togetehr two grids with values of "
00459                         "different multiplicity.\n");
00460     if (scale_factor != 1.0) 
00461       for (size_t i = 0; i < data.size(); i++) {
00462         data[i] += scale_factor * other_grid.data[i];
00463       }
00464     else 
00465       // skip multiplication if possible
00466       for (size_t i = 0; i < data.size(); i++) {
00467         data[i] += other_grid.data[i];
00468       }
00469     has_data = true;
00470   }
00471 
00474   virtual inline T value_output (std::vector<int> const &ix,
00475                                  size_t const &imult = 0)
00476   {
00477     return value (ix, imult);
00478   }
00479 
00483   virtual inline void value_input (std::vector<int> const &ix,
00484                                    T const &t,
00485                                    size_t const &imult = 0,
00486                                    bool add = false)
00487   {
00488     if ( add )
00489       data[address (ix) + imult] += t;
00490     else
00491       data[address (ix) + imult] = t;
00492     has_data = true;
00493   }
00494 
00495   //   /// Get the pointer to the binned value indexed by ix
00496   //   inline T const *value_p (std::vector<int> const &ix)
00497   //   {
00498   //     return &(data[address (ix)]);
00499   //   }
00500 
00503   inline std::vector<int> const new_index() const
00504   {
00505     return std::vector<int> (nd, 0);
00506   }
00507 
00510   inline bool index_ok (std::vector<int> const &ix) const
00511   {
00512     for (size_t i = 0; i < nd; i++) {
00513       if ( (ix[i] < 0) || (ix[i] >= int (nx[i])) )
00514         return false;
00515     }
00516     return true;
00517   }
00518 
00521   inline void incr (std::vector<int> &ix) const
00522   {
00523     for (int i = ix.size()-1; i >= 0; i--) {
00524 
00525       ix[i]++;
00526 
00527       if (ix[i] >= nx[i]) {
00528 
00529         if (i > 0) {
00530           ix[i] = 0;
00531           continue;
00532         } else {
00533           // this is the last iteration, a non-valid index is being
00534           // set for the outer index, which will be caught by
00535           // index_ok()
00536           ix[0] = nx[0];
00537           return;
00538         }
00539       } else {
00540         return;
00541       }
00542     }
00543   }
00544 
00546   std::ostream & write_params (std::ostream &os)
00547   {
00548     os << "grid_parameters {\n  n_colvars " << nd << "\n";
00549       
00550     os << "  lower_boundaries ";
00551     for (size_t i = 0; i < nd; i++) 
00552       os << " " << lower_boundaries[i];
00553     os << "\n";
00554 
00555     os << "  upper_boundaries ";
00556     for (size_t i = 0; i < nd; i++) 
00557       os << " " << upper_boundaries[i];
00558     os << "\n";
00559 
00560     os << "  widths ";
00561     for (size_t i = 0; i < nd; i++) 
00562       os << " " << widths[i];
00563     os << "\n";
00564 
00565     os << "  sizes ";
00566     for (size_t i = 0; i < nd; i++) 
00567       os << " " << nx[i];
00568     os << "\n";
00569 
00570     os << "}\n";
00571     return os;
00572   } 
00573 
00574 
00575   bool parse_params (std::string const &conf)
00576   {
00577     std::vector<int> old_nx = nx;
00578     std::vector<colvarvalue> old_lb = lower_boundaries;
00579 
00580     {
00581       size_t nd_in = 0;
00582       colvarparse::get_keyval (conf, "n_colvars", nd_in, nd, colvarparse::parse_silent);
00583       if (nd_in != nd)
00584         cvm::fatal_error ("Error: trying to read data for a grid "
00585                           "that contains a different number of colvars ("+
00586                           cvm::to_str (nd_in)+") than the grid defined "
00587                           "in the configuration file ("+cvm::to_str (nd)+
00588                           ").\n");
00589     }
00590 
00591     colvarparse::get_keyval (conf, "lower_boundaries",
00592                              lower_boundaries, lower_boundaries, colvarparse::parse_silent);
00593 
00594     colvarparse::get_keyval (conf, "upper_boundaries",
00595                              upper_boundaries, upper_boundaries, colvarparse::parse_silent);
00596 
00597     colvarparse::get_keyval (conf, "widths", widths, widths, colvarparse::parse_silent);
00598 
00599     colvarparse::get_keyval (conf, "sizes", nx, nx, colvarparse::parse_silent);
00600 
00601     bool new_params = false;
00602     for (size_t i = 0; i < nd; i++) {
00603       if ( (old_nx[i] != nx[i]) ||
00604            (std::sqrt (cv[i]->dist2 (old_lb[i],
00605                                      lower_boundaries[i])) > 1.0E-10) ) {
00606         new_params = true;
00607       }
00608     }
00609 
00610     // reallocate the array in case the grid params have just changed
00611     if (new_params) {
00612       data.resize (0);
00613       this->create (nx, T(), mult);
00614     }
00615 
00616     return true;
00617   }
00618 
00622   void check_consistency()
00623   {
00624     for (size_t i = 0; i < nd; i++) {
00625       if ( (std::sqrt (cv[i]->dist2 (cv[i]->lower_boundary,
00626                                      lower_boundaries[i])) > 1.0E-10) || 
00627            (std::sqrt (cv[i]->dist2 (cv[i]->upper_boundary,
00628                                      upper_boundaries[i])) > 1.0E-10) || 
00629            (std::sqrt (cv[i]->dist2 (cv[i]->width,
00630                                      widths[i])) > 1.0E-10) ) {
00631         cvm::fatal_error ("Error: restart information for a grid is "
00632                           "inconsistent with that of its colvars.\n");
00633       }
00634     }
00635   }
00636 
00637 
00640   void check_consistency (colvar_grid<T> const &other_grid)
00641   {
00642     for (size_t i = 0; i < nd; i++) {
00643       // we skip dist2(), because periodicities and the like should
00644       // matter: boundaries should be EXACTLY the same (otherwise,
00645       // map_grid() should be used)
00646       if ( (std::fabs (other_grid.lower_boundaries[i] -
00647                        lower_boundaries[i]) > 1.0E-10) || 
00648            (std::fabs (other_grid.upper_boundaries[i] -
00649                        upper_boundaries[i]) > 1.0E-10) || 
00650            (std::fabs (other_grid.widths[i] -
00651                        widths[i]) > 1.0E-10) || 
00652            (data.size() != other_grid.data.size()) ) {
00653       cvm::fatal_error ("Error: inconsistency between "
00654                         "two grids that are supposed to be equal, "
00655                         "aside from the data stored.\n");
00656     }
00657   }
00658 }
00659   
00660 
00664   std::ostream & write_raw (std::ostream &os,
00665                             size_t const buf_size = 3)
00666   {
00667     std::streamsize const w = os.width();
00668     std::streamsize const p = os.precision();
00669 
00670     std::vector<int> ix = new_index();
00671     size_t count = 0;
00672     for ( ; index_ok (ix); incr (ix)) {
00673       for (size_t imult = 0; imult < mult; imult++) {
00674         os << " "
00675            << std::setw (w) << std::setprecision (p)
00676            << value_output (ix, imult);
00677         if (((++count) % buf_size) == 0)
00678           os << "\n";
00679       }
00680     }
00681     // write a final newline only if buffer is not empty
00682     if ((count % buf_size) != 0)
00683       os << "\n";
00684 
00685     return os;
00686   }
00687 
00689 std::istream & read_raw (std::istream &is)
00690 {
00691   size_t const start_pos = is.tellg();
00692     
00693   for (std::vector<int> ix = new_index(); index_ok (ix); incr (ix)) {
00694     for (size_t imult = 0; imult < mult; imult++) {
00695       T new_value;
00696       if (is >> new_value) {
00697         value_input (ix, new_value, imult);
00698       } else {
00699         is.clear();
00700         is.seekg (start_pos, std::ios::beg);
00701         is.setstate (std::ios::failbit);
00702         return is;
00703       }
00704     }
00705   }
00706 
00707   has_data = true;
00708   return is;
00709 }
00710 
00712 void read_raw_error()
00713 {
00714   cvm::fatal_error ("Error: failed to read all of the grid points from file.  Possible explanations: grid parameters in the configuration (lowerBoundary, upperBoundary, width) are different from those in the file, or the file is corrupt/incomplete.\n");
00715 }
00716 
00719 void write_multicol (std::ostream &os)
00720 {
00721   std::streamsize const w = os.width();
00722   std::streamsize const p = os.precision();
00723 
00724   // Data in the header: nColvars, then for each
00725   // xiMin, dXi, nPoints, periodic
00726 
00727   os << std::setw (2) << "# " << nd << "\n";
00728   for (size_t i = 0; i < nd; i++) {
00729     os << "# "
00730        << std::setw (10) << lower_boundaries[i]
00731        << std::setw (10) << widths[i]
00732        << std::setw (10) << nx[i] << "  "
00733        << periodic[i] << "\n";
00734   }
00735 
00736   for (std::vector<int> ix = new_index(); index_ok (ix); incr (ix) ) {
00737 
00738     if (ix.back() == 0) {
00739       // if the last index is 0, add a new line to mark the new record
00740       os << "\n";
00741     }
00742 
00743     for (size_t i = 0; i < nd; i++) {
00744       os << " "
00745          << std::setw (w) << std::setprecision (p)
00746          << bin_to_value_scalar (ix[i], i);
00747     }
00748     os << " ";
00749     for (size_t imult = 0; imult < mult; imult++) {
00750       os << " "
00751          << std::setw (w) << std::setprecision (p)
00752          << value_output (ix, imult);
00753     }
00754     os << "\n";
00755   }
00756 }
00757 
00760 std::istream & read_multicol (std::istream &is, bool add = false)
00761 {
00762   // Data in the header: nColvars, then for each
00763   // xiMin, dXi, nPoints, periodic
00764 
00765   std::string   hash;
00766   cvm::real     lower, width, x;
00767   size_t        n, periodic;
00768   bool          remap;
00769   std::vector<T>        new_value;
00770   std::vector<int>      nx_read;
00771   std::vector<int>      bin;
00772 
00773   if ( cv.size() != nd ) {
00774     cvm::fatal_error ("Cannot read grid file: missing reference to colvars.");
00775   }
00776 
00777   if ( !(is >> hash) || (hash != "#") ) {
00778     cvm::fatal_error ("Error reading grid at position "+
00779                       cvm::to_str (is.tellg())+" in stream (read \"" + hash + "\")\n");
00780   }
00781 
00782   is >> n;
00783   if ( n != nd ) {
00784     cvm::fatal_error ("Error reading grid: wrong number of collective variables.\n");
00785   }
00786 
00787   nx_read.resize (n);
00788   bin.resize (n);
00789   new_value.resize (mult);
00790 
00791   if (this->has_parent_data && add) {
00792     new_data.resize (data.size());
00793   }
00794 
00795   remap = false;
00796   for (size_t i = 0; i < nd; i++ ) {
00797     if ( !(is >> hash) || (hash != "#") ) {
00798       cvm::fatal_error ("Error reading grid at position "+
00799                         cvm::to_str (is.tellg())+" in stream (read \"" + hash + "\")\n");
00800     }
00801 
00802     is >> lower >> width >> nx_read[i] >> periodic;
00803 
00804 
00805     if ( (std::fabs (lower - lower_boundaries[i].real_value) > 1.0e-10) ||
00806          (std::fabs (width - widths[i] ) > 1.0e-10) ||
00807          (nx_read[i] != nx[i]) ) {
00808       cvm::log ("Warning: reading from different grid definition (colvar "
00809                 + cvm::to_str (i+1) + "); remapping data on new grid.\n");
00810       remap = true;
00811     }
00812   }
00813 
00814   if ( remap ) {
00815     // re-grid data
00816     while (is.good()) {
00817       bool end_of_file = false;
00818 
00819       for (size_t i = 0; i < nd; i++ ) {
00820         if ( !(is >> x) ) end_of_file = true;
00821         bin[i] = value_to_bin_scalar (x, i);
00822       }
00823       if (end_of_file) break;   
00824 
00825       for (size_t imult = 0; imult < mult; imult++) {
00826         is >> new_value[imult];
00827       }
00828 
00829       if ( index_ok(bin) ) {
00830         for (size_t imult = 0; imult < mult; imult++) {
00831           value_input (bin, new_value[imult], imult, add);
00832         }
00833       }
00834     }
00835   } else {
00836     // do not re-grid the data but assume the same grid is used
00837     for (std::vector<int> ix = new_index(); index_ok (ix); incr (ix) ) {
00838       for (size_t i = 0; i < nd; i++ ) {
00839         is >> x;
00840       }
00841       for (size_t imult = 0; imult < mult; imult++) {
00842         is >> new_value[imult];
00843         value_input (ix, new_value[imult], imult, add);
00844       }
00845     }
00846   }
00847   has_data = true;
00848   return is;
00849 }
00850 
00851 };
00852 
00853 
00854 
00857 class colvar_grid_count : public colvar_grid<size_t>
00858 {
00859 public:
00860 
00862   colvar_grid_count();
00863 
00865   virtual inline ~colvar_grid_count()
00866   {}
00867 
00869   colvar_grid_count (std::vector<int> const &nx_i,
00870                      size_t const           &def_count = 0);
00871 
00873   colvar_grid_count (std::vector<colvar *>  &colvars,
00874                      size_t const           &def_count = 0);
00875 
00877   inline void incr_count (std::vector<int> const &ix)
00878   {
00879     ++(data[this->address (ix)]);
00880   }
00881 
00883   inline size_t const & new_count (std::vector<int> const &ix,
00884                                    size_t const &imult = 0)
00885   {
00886     return new_data[address (ix) + imult];
00887   }
00888 
00890   std::istream & read_restart (std::istream &is);
00891 
00893   std::ostream & write_restart (std::ostream &os);
00894 
00898   virtual inline void value_input (std::vector<int> const &ix,
00899                                    size_t const &t,
00900                                    size_t const &imult = 0,
00901                                    bool add = false)
00902   {
00903     if (add) {
00904       data[address (ix)] += t;
00905       if (this->has_parent_data) {
00906         // save newly read data for inputting parent grid
00907         new_data[address (ix)] = t;
00908       }
00909     } else {
00910       data[address (ix)] = t;
00911     }
00912     has_data = true;
00913   }
00914 };
00915 
00916 
00918 class colvar_grid_scalar : public colvar_grid<cvm::real>
00919 {
00920 public:
00921 
00924   colvar_grid_count *samples;
00925 
00927   colvar_grid_scalar();
00928 
00930   colvar_grid_scalar (colvar_grid_scalar const &g);
00931 
00933   ~colvar_grid_scalar();
00934 
00936   colvar_grid_scalar (std::vector<int> const &nx_i);
00937 
00939   colvar_grid_scalar (std::vector<colvar *> &colvars,
00940                       bool margin = 0);
00941 
00943   inline void acc_value (std::vector<int> const &ix,
00944                          cvm::real const &new_value,
00945                          size_t const &imult = 0)
00946   {
00947     // only legal value of imult here is 0
00948     data[address (ix)] += new_value;
00949     if (samples)
00950       samples->incr_count (ix);
00951     has_data = true;
00952   }
00953 
00955   inline const cvm::real * gradient_finite_diff ( const std::vector<int> &ix0 )
00956   {
00957     cvm::real A0, A1;
00958     std::vector<int> ix;
00959     if (nd != 2) cvm::fatal_error ("Finite differences available in dimension 2 only.");
00960     for (int n = 0; n < nd; n++) {
00961       ix = ix0;
00962       A0 = data[address (ix)];
00963       ix[n]++; wrap (ix);
00964       A1 = data[address (ix)];
00965       ix[1-n]++; wrap (ix);
00966       A1 += data[address (ix)];
00967       ix[n]--; wrap (ix);
00968       A0 += data[address (ix)];
00969       grad[n] = 0.5 * (A1 - A0) / widths[n];
00970     }
00971     return grad; 
00972   }
00973 
00976   virtual cvm::real value_output (std::vector<int> const &ix,
00977                                   size_t const &imult = 0)
00978   {
00979     if (imult > 0)
00980       cvm::fatal_error ("Error: trying to access a component "
00981                         "larger than 1 in a scalar data grid.\n");
00982     if (samples)
00983       return (samples->value (ix) > 0) ?
00984         (data[address (ix)] / cvm::real (samples->value (ix))) :
00985         0.0;
00986     else
00987       return data[address (ix)];
00988   }
00989 
00993   virtual void value_input (std::vector<int> const &ix,
00994                             cvm::real const &new_value,
00995                             size_t const &imult = 0,
00996                             bool add = false)
00997   {
00998     if (imult > 0)
00999       cvm::fatal_error ("Error: trying to access a component "
01000                         "larger than 1 in a scalar data grid.\n");
01001     if (add) {
01002       if (samples)
01003         data[address (ix)] += new_value * samples->new_count (ix);
01004       else
01005         data[address (ix)] += new_value;
01006     } else {
01007       if (samples)
01008         data[address (ix)] = new_value * samples->value (ix);
01009       else
01010         data[address (ix)] = new_value;
01011     }
01012     has_data = true;
01013   }
01014 
01016   std::istream & read_restart (std::istream &is);
01017 
01019   std::ostream & write_restart (std::ostream &os);
01020 
01022   inline cvm::real maximum_value()
01023   {
01024     cvm::real max = data[0];
01025     for (size_t i = 0; i < nt; i++) {
01026       if (data[i] > max) max = data[i];
01027     }
01028     return max;
01029   }
01030 
01032   inline cvm::real minimum_value()
01033   {
01034     cvm::real min = data[0];
01035     for (size_t i = 0; i < nt; i++) {
01036       if (data[i] < min) min = data[i];
01037     }
01038     return min;
01039   }
01040 
01041 private:
01042   // gradient
01043   cvm::real * grad;
01044 };
01045 
01046 
01047 
01049 class colvar_grid_gradient : public colvar_grid<cvm::real>
01050 {
01051 public:
01052 
01055   colvar_grid_count *samples;
01056 
01058   colvar_grid_gradient();
01059 
01061   virtual inline ~colvar_grid_gradient()
01062   {}
01063 
01065   colvar_grid_gradient (std::vector<int> const &nx_i);
01066 
01068   colvar_grid_gradient (std::vector<colvar *>  &colvars);
01069 
01071   inline void acc_grad (std::vector<int> const &ix, cvm::real const *grads) {
01072     for (size_t imult = 0; imult < mult; imult++) {
01073       data[address (ix) + imult] += grads[imult];
01074     }
01075     if (samples)
01076       samples->incr_count (ix);
01077   }
01078 
01081   inline void acc_force (std::vector<int> const &ix, cvm::real const *forces) {
01082     for (size_t imult = 0; imult < mult; imult++) {
01083       data[address (ix) + imult] -= forces[imult];
01084     }
01085     if (samples)
01086       samples->incr_count (ix);
01087   }
01088 
01091   virtual inline cvm::real value_output (std::vector<int> const &ix,
01092                                          size_t const &imult = 0)
01093   {
01094     if (samples)
01095       return (samples->value (ix) > 0) ?
01096         (data[address (ix) + imult] / cvm::real (samples->value (ix))) :
01097         0.0;
01098     else
01099       return data[address (ix) + imult];
01100   }
01101 
01105   virtual inline void value_input (std::vector<int> const &ix,
01106                                    cvm::real const &new_value,
01107                                    size_t const &imult = 0,
01108                                    bool add = false)
01109   {
01110     if (add) {
01111       if (samples)
01112         data[address (ix) + imult] += new_value * samples->new_count (ix);
01113       else
01114         data[address (ix) + imult] += new_value;
01115     } else {
01116       if (samples)
01117         data[address (ix) + imult] = new_value * samples->value (ix);
01118       else
01119         data[address (ix) + imult] = new_value;
01120     }
01121     has_data = true;
01122   }
01123 
01124 
01126   std::istream & read_restart (std::istream &is);
01127 
01129   std::ostream & write_restart (std::ostream &os);
01130 
01132   inline cvm::real average()
01133   {
01134     size_t n = 0;
01135 
01136     if (nd != 1 || nx[0] == 0) {
01137       return 0.0;
01138     }
01139 
01140     cvm::real sum = 0.0;
01141     std::vector<int> ix = new_index();
01142     if (samples) {
01143       for ( ; index_ok (ix); incr (ix)) {
01144         if ( (n = samples->value (ix)) )
01145           sum += value (ix) / n;
01146       }
01147     } else {
01148       for ( ; index_ok (ix); incr (ix)) {
01149         sum += value (ix);
01150       }
01151     }
01152     return (sum / cvm::real (nx[0]));
01153   }
01154 
01157   void write_1D_integral (std::ostream &os);
01158 
01159 };
01160 
01161 
01162 #endif
01163 

Generated on Fri May 25 04:07:13 2012 for NAMD by  doxygen 1.3.9.1