Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2016-2021 The VES code team
       3             :    (see the PEOPLE-VES file at the root of this folder for a list of names)
       4             : 
       5             :    See http://www.ves-code.org for more information.
       6             : 
       7             :    This file is part of VES code module.
       8             : 
       9             :    The VES code module is free software: you can redistribute it and/or modify
      10             :    it under the terms of the GNU Lesser General Public License as published by
      11             :    the Free Software Foundation, either version 3 of the License, or
      12             :    (at your option) any later version.
      13             : 
      14             :    The VES code module is distributed in the hope that it will be useful,
      15             :    but WITHOUT ANY WARRANTY; without even the implied warranty of
      16             :    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
      17             :    GNU Lesser General Public License for more details.
      18             : 
      19             :    You should have received a copy of the GNU Lesser General Public License
      20             :    along with the VES code module.  If not, see <http://www.gnu.org/licenses/>.
      21             : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
      22             : 
      23             : #include "CoeffsMatrix.h"
      24             : #include "CoeffsVector.h"
      25             : #include "BasisFunctions.h"
      26             : 
      27             : #include "tools/Tools.h"
      28             : #include "core/Value.h"
      29             : #include "tools/File.h"
      30             : #include "tools/Exception.h"
      31             : #include "tools/Random.h"
      32             : #include "tools/Communicator.h"
      33             : 
      34             : #include <vector>
      35             : #include <cmath>
      36             : #include <iostream>
      37             : #include <sstream>
      38             : #include <cstdio>
      39             : #include <cfloat>
      40             : 
      41             : 
      42             : namespace PLMD {
      43             : namespace ves {
      44             : 
      45           0 : CoeffsMatrix::CoeffsMatrix(
      46             :   const std::string& label,
      47             :   const std::vector<std::string>& dimension_labels,
      48             :   const std::vector<unsigned int>& indices_shape,
      49             :   Communicator& cc,
      50             :   const bool diagonal,
      51           0 :   const bool use_iteration_counter):
      52             :   CoeffsBase(label,dimension_labels,indices_shape,use_iteration_counter),
      53           0 :   data(0),
      54           0 :   size_(0),
      55           0 :   nrows_(0),
      56           0 :   ncolumns_(0),
      57           0 :   diagonal_(diagonal),
      58           0 :   averaging_counter(0),
      59           0 :   averaging_exp_decay_(0),
      60           0 :   mycomm(cc)
      61             : {
      62           0 :   setupMatrix();
      63           0 : }
      64             : 
      65             : 
      66           0 : CoeffsMatrix::CoeffsMatrix(
      67             :   const std::string& label,
      68             :   std::vector<Value*>& args,
      69             :   std::vector<BasisFunctions*>& basisf,
      70             :   Communicator& cc,
      71             :   const bool diagonal,
      72           0 :   const bool use_iteration_counter):
      73             :   CoeffsBase(label,args,basisf,use_iteration_counter),
      74           0 :   data(0),
      75           0 :   size_(0),
      76           0 :   nrows_(0),
      77           0 :   ncolumns_(0),
      78           0 :   diagonal_(diagonal),
      79           0 :   averaging_counter(0),
      80           0 :   averaging_exp_decay_(0),
      81           0 :   mycomm(cc)
      82             : {
      83           0 :   setupMatrix();
      84           0 : }
      85             : 
      86             : 
      87           0 : CoeffsMatrix::CoeffsMatrix(
      88             :   const std::string& label,
      89             :   std::vector<std::vector<Value*> >& argsv,
      90             :   std::vector<std::vector<BasisFunctions*> >& basisfv,
      91             :   Communicator& cc,
      92             :   const bool diagonal,
      93             :   const bool use_iteration_counter,
      94           0 :   const std::string& multicoeffs_label):
      95             :   CoeffsBase(label,argsv,basisfv,use_iteration_counter,multicoeffs_label),
      96           0 :   data(0),
      97           0 :   size_(0),
      98           0 :   nrows_(0),
      99           0 :   ncolumns_(0),
     100           0 :   diagonal_(diagonal),
     101           0 :   averaging_counter(0),
     102           0 :   averaging_exp_decay_(0),
     103           0 :   mycomm(cc)
     104             : {
     105           0 :   setupMatrix();
     106           0 : }
     107             : 
     108             : 
     109         172 : CoeffsMatrix::CoeffsMatrix(
     110             :   const std::string& label,
     111             :   CoeffsVector* coeffsVec,
     112             :   Communicator& cc,
     113         172 :   const bool diagonal):
     114             :   CoeffsBase( *(static_cast<CoeffsBase*>(coeffsVec)) ),
     115         172 :   data(0),
     116         172 :   size_(0),
     117         172 :   nrows_(0),
     118         172 :   ncolumns_(0),
     119         172 :   diagonal_(diagonal),
     120         172 :   averaging_counter(0),
     121         172 :   averaging_exp_decay_(0),
     122         172 :   mycomm(cc)
     123             : {
     124         172 :   setLabels(label);
     125         172 :   setupMatrix();
     126         172 : }
     127             : 
     128             : 
     129         172 : CoeffsMatrix::~CoeffsMatrix() {}
     130             : 
     131             : 
     132         172 : void CoeffsMatrix::setupMatrix() {
     133         172 :   nrows_=numberOfCoeffs();
     134         172 :   ncolumns_=nrows_;
     135         172 :   if(diagonal_) {
     136         172 :     size_=nrows_;
     137             :   }
     138             :   else {
     139           0 :     size_=(nrows_*nrows_-nrows_)/2+nrows_;
     140             :   }
     141         172 :   clear();
     142         172 : }
     143             : 
     144             : 
     145         439 : size_t CoeffsMatrix::getSize() const {
     146         439 :   return size_;
     147             : }
     148             : 
     149             : 
     150       23395 : bool CoeffsMatrix::isDiagonal() const {
     151       23395 :   return diagonal_;
     152             : }
     153             : 
     154             : 
     155           0 : bool CoeffsMatrix::sameShape(CoeffsVector& coeffsvector_in) const {
     156           0 :   return CoeffsBase::sameShape( *(static_cast<CoeffsBase*>(&coeffsvector_in)) );
     157             : }
     158             : 
     159             : 
     160           0 : bool CoeffsMatrix::sameShape(CoeffsMatrix& coeffsmat_in) const {
     161           0 :   return CoeffsBase::sameShape( *(static_cast<CoeffsBase*>(&coeffsmat_in)) );
     162             : }
     163             : 
     164             : 
     165           0 : bool CoeffsMatrix::sameShape(CoeffsMatrix& coeffsmat0, CoeffsMatrix& coeffsmat1) {
     166           0 :   return coeffsmat0.sameShape(coeffsmat1);
     167             : }
     168             : 
     169             : 
     170           0 : bool CoeffsMatrix::sameShape(CoeffsVector& coeffsvec, CoeffsMatrix& coeffsmat) {
     171           0 :   return coeffsmat.sameShape(coeffsvec);
     172             : }
     173             : 
     174             : 
     175           0 : bool CoeffsMatrix::sameShape(CoeffsMatrix& coeffsmat, CoeffsVector& coeffsvec) {
     176           0 :   return coeffsmat.sameShape(coeffsvec);
     177             : }
     178             : 
     179             : 
     180           0 : void CoeffsMatrix::sumCommMPI() {
     181           0 :   mycomm.Sum(data);
     182           0 : }
     183             : 
     184             : 
     185           0 : void CoeffsMatrix::sumCommMPI(Communicator& cc) {
     186           0 :   cc.Sum(data);
     187           0 : }
     188             : 
     189             : 
     190           0 : void CoeffsMatrix::sumMultiSimCommMPI(Communicator& multi_sim_cc) {
     191           0 :   if(mycomm.Get_rank()==0) {
     192           0 :     double nwalkers = static_cast<double>(multi_sim_cc.Get_size());
     193           0 :     multi_sim_cc.Sum(data);
     194           0 :     scaleAllValues(1.0/nwalkers);
     195             :   }
     196           0 :   mycomm.Bcast(data,0);
     197           0 : }
     198             : 
     199             : 
     200     5425906 : size_t CoeffsMatrix::getMatrixIndex(const size_t index1, const size_t index2) const {
     201             :   size_t matrix_idx;
     202             :   plumed_dbg_assert(index1<nrows_);
     203             :   plumed_dbg_assert(index2<ncolumns_);
     204     5425906 :   if(diagonal_) {
     205             :     // plumed_massert(index1==index2,"CoeffsMatrix: you trying to access a off-diagonal element of a diagonal coeffs matrix");
     206             :     matrix_idx=index1;
     207             :   }
     208           0 :   else if (index1<=index2) {
     209           0 :     matrix_idx=index2+index1*(nrows_-1)-index1*(index1-1)/2;
     210             :   }
     211             :   else {
     212           0 :     matrix_idx=index1+index2*(nrows_-1)-index2*(index2-1)/2;
     213             :   }
     214     5425906 :   return matrix_idx;
     215             : }
     216             : 
     217             : 
     218         267 : void CoeffsMatrix::clear() {
     219         267 :   data.resize(getSize());
     220       19388 :   for(size_t i=0; i<data.size(); i++) {
     221       19121 :     data[i]=0.0;
     222             :   }
     223         267 : }
     224             : 
     225             : 
     226           0 : void CoeffsMatrix::setAllValuesToZero() {
     227           0 :   for(size_t i=0; i<data.size(); i++) {
     228           0 :     data[i]=0.0;
     229             :   }
     230           0 : }
     231             : 
     232             : 
     233       27630 : double CoeffsMatrix::getValue(const size_t index1, const size_t index2) const {
     234       27630 :   return data[getMatrixIndex(index1,index2)];
     235             : }
     236             : 
     237             : 
     238           0 : double CoeffsMatrix::getValue(const std::vector<unsigned int>& indices1, const std::vector<unsigned int>& indices2) const {
     239           0 :   return getValue(getIndex(indices1),getIndex(indices2));
     240             : }
     241             : 
     242             : 
     243         219 : void CoeffsMatrix::setValue(const size_t index1, const size_t index2, const double value) {
     244         219 :   data[getMatrixIndex(index1,index2)]=value;
     245         219 : }
     246             : 
     247             : 
     248           0 : void CoeffsMatrix::setValue(const std::vector<unsigned int>& indices1, const std::vector<unsigned int>& indices2, const double value) {
     249           0 :   setValue(getIndex(indices1),getIndex(indices2),value);
     250           0 : }
     251             : 
     252             : 
     253           0 : double& CoeffsMatrix::operator()(const size_t index1, const size_t index2) {
     254           0 :   return data[getMatrixIndex(index1,index2)];
     255             : }
     256             : 
     257             : 
     258     1790345 : const double& CoeffsMatrix::operator()(const size_t index1, const size_t index2) const {
     259     1790345 :   return data[getMatrixIndex(index1,index2)];
     260             : }
     261             : 
     262             : 
     263           0 : double& CoeffsMatrix::operator()(const std::vector<unsigned int>& indices1, const std::vector<unsigned int>& indices2) {
     264           0 :   return data[getMatrixIndex(getIndex(indices1),getIndex(indices2))];
     265             : }
     266             : 
     267             : 
     268           0 : const double& CoeffsMatrix::operator()(const std::vector<unsigned int>& indices1, const std::vector<unsigned int>& indices2) const {
     269           0 :   return data[getMatrixIndex(getIndex(indices1),getIndex(indices2))];
     270             : }
     271             : 
     272             : 
     273       22735 : CoeffsVector operator*(const CoeffsMatrix& coeffs_matrix, const CoeffsVector& coeffs_vector) {
     274       22735 :   CoeffsVector new_coeffs_vector(coeffs_vector);
     275       22735 :   new_coeffs_vector.clear();
     276       22735 :   plumed_massert(coeffs_vector.numberOfCoeffs()==coeffs_matrix.numberOfCoeffs(),"CoeffsMatrix and CoeffsVector are of the wrong size");
     277             :   size_t numcoeffs = coeffs_vector.numberOfCoeffs();
     278       22735 :   if(coeffs_matrix.isDiagonal()) {
     279     1813080 :     for(size_t i=0; i<numcoeffs; i++) {
     280     1790345 :       new_coeffs_vector(i) = coeffs_matrix(i,i)*coeffs_vector(i);
     281             :     }
     282             :   }
     283             :   else {
     284           0 :     for(size_t i=0; i<numcoeffs; i++) {
     285           0 :       for(size_t j=0; j<numcoeffs; j++) {
     286           0 :         new_coeffs_vector(i) += coeffs_matrix(i,j)*coeffs_vector(j);
     287             :       }
     288             :     }
     289             :   }
     290       22735 :   return new_coeffs_vector;
     291           0 : }
     292             : 
     293             : 
     294           0 : void CoeffsMatrix::addToValue(const size_t index1, const size_t index2, const double value) {
     295           0 :   data[getMatrixIndex(index1,index2)]+=value;
     296           0 : }
     297             : 
     298             : 
     299           0 : void CoeffsMatrix::addToValue(const std::vector<unsigned int>& indices1, const std::vector<unsigned int>& indices2, const double value) {
     300           0 :   addToValue(getIndex(indices1),getIndex(indices2),value);
     301           0 : }
     302             : 
     303             : 
     304       22810 : void CoeffsMatrix::scaleAllValues(const double scalef) {
     305     1823850 :   for(size_t i=0; i<data.size(); i++) {
     306     1801040 :     data[i]*=scalef;
     307             :   }
     308       22810 : }
     309             : 
     310             : 
     311       22810 : CoeffsMatrix& CoeffsMatrix::operator*=(const double scalef) {
     312       22810 :   scaleAllValues(scalef);
     313       22810 :   return *this;
     314             : }
     315             : 
     316             : 
     317           0 : CoeffsMatrix operator*(const double scalef, const CoeffsMatrix& coeffsmatrix) {
     318           0 :   return CoeffsMatrix(coeffsmatrix)*=scalef;
     319             : }
     320             : 
     321             : 
     322           0 : CoeffsMatrix operator*(const CoeffsMatrix& coeffsmatrix, const double scalef) {
     323           0 :   return scalef*coeffsmatrix;
     324             : }
     325             : 
     326             : 
     327           0 : CoeffsMatrix& CoeffsMatrix::operator*=(const CoeffsMatrix& other_coeffsmatrix) {
     328           0 :   plumed_massert(data.size()==other_coeffsmatrix.getSize(),"Coeffs matrices do not have the same size");
     329           0 :   for(size_t i=0; i<data.size(); i++) {
     330           0 :     data[i]*=other_coeffsmatrix.data[i];
     331             :   }
     332           0 :   return *this;
     333             : }
     334             : 
     335             : 
     336           0 : CoeffsMatrix CoeffsMatrix::operator*(const CoeffsMatrix& other_coeffsmatrix) const {
     337           0 :   return CoeffsMatrix(*this)*=other_coeffsmatrix;
     338             : }
     339             : 
     340             : 
     341           0 : void CoeffsMatrix::setValues(const double value) {
     342           0 :   for(size_t i=0; i<data.size(); i++) {
     343           0 :     data[i]=value;
     344             :   }
     345           0 : }
     346             : 
     347             : 
     348       22810 : void CoeffsMatrix::setValues(const std::vector<double>& values) {
     349       22810 :   plumed_massert( data.size()==values.size(), "Incorrect size");
     350     1823850 :   for(size_t i=0; i<data.size(); i++) {
     351     1801040 :     data[i]=values[i];
     352             :   }
     353       22810 : }
     354             : 
     355             : 
     356           0 : void CoeffsMatrix::setValues(const CoeffsMatrix& other_coeffsmatrix) {
     357           0 :   plumed_massert( data.size()==other_coeffsmatrix.getSize(), "Incorrect size");
     358           0 :   for(size_t i=0; i<data.size(); i++) {
     359           0 :     data[i]=other_coeffsmatrix.data[i];
     360             :   }
     361           0 : }
     362             : 
     363             : 
     364           0 : CoeffsMatrix& CoeffsMatrix::operator=(const double value) {
     365           0 :   setValues(value);
     366           0 :   return *this;
     367             : }
     368             : 
     369             : 
     370       22810 : CoeffsMatrix& CoeffsMatrix::operator=(const std::vector<double>& values) {
     371       22810 :   setValues(values);
     372       22810 :   return *this;
     373             : }
     374             : 
     375             : 
     376             : // CoeffsMatrix& CoeffsMatrix::operator=(const CoeffsMatrix& other_coeffsmatrix) {
     377             : //   setValues(other_coeffsmatrix);
     378             : //   return *this;
     379             : // }
     380             : 
     381             : 
     382           0 : CoeffsMatrix CoeffsMatrix::operator+() const {
     383           0 :   return *this;
     384             : }
     385             : 
     386             : 
     387           0 : CoeffsMatrix CoeffsMatrix::operator-() const {
     388           0 :   return CoeffsMatrix(*this)*=-1.0;
     389             : }
     390             : 
     391             : 
     392           0 : void CoeffsMatrix::addToValues(const double value) {
     393           0 :   for(size_t i=0; i<data.size(); i++) {
     394           0 :     data[i]+=value;
     395             :   }
     396           0 : }
     397             : 
     398             : 
     399           0 : void CoeffsMatrix::addToValues(const std::vector<double>& values) {
     400           0 :   plumed_massert( data.size()==values.size(), "Incorrect size");
     401           0 :   for(size_t i=0; i<data.size(); i++) {
     402           0 :     data[i]+=values[i];
     403             :   }
     404           0 : }
     405             : 
     406             : 
     407           0 : void CoeffsMatrix::addToValues(const CoeffsMatrix& other_coeffsmatrix) {
     408           0 :   plumed_massert( data.size()==other_coeffsmatrix.getSize(), "Incorrect size");
     409           0 :   for(size_t i=0; i<data.size(); i++) {
     410           0 :     data[i]+=other_coeffsmatrix.data[i];
     411             :   }
     412           0 : }
     413             : 
     414             : 
     415           0 : void CoeffsMatrix::subtractFromValues(const double value) {
     416           0 :   for(size_t i=0; i<data.size(); i++) {
     417           0 :     data[i]-=value;
     418             :   }
     419           0 : }
     420             : 
     421             : 
     422           0 : void CoeffsMatrix::subtractFromValues(const std::vector<double>& values) {
     423           0 :   plumed_massert( data.size()==values.size(), "Incorrect size");
     424           0 :   for(size_t i=0; i<data.size(); i++) {
     425           0 :     data[i]-=values[i];
     426             :   }
     427           0 : }
     428             : 
     429             : 
     430           0 : void CoeffsMatrix::subtractFromValues(const CoeffsMatrix& other_coeffsmatrix) {
     431           0 :   plumed_massert( data.size()==other_coeffsmatrix.getSize(), "Incorrect size");
     432           0 :   for(size_t i=0; i<data.size(); i++) {
     433           0 :     data[i]-=other_coeffsmatrix.data[i];
     434             :   }
     435           0 : }
     436             : 
     437             : 
     438           0 : CoeffsMatrix& CoeffsMatrix::operator+=(const double value) {
     439           0 :   addToValues(value);
     440           0 :   return *this;
     441             : }
     442             : 
     443             : 
     444           0 : CoeffsMatrix operator+(const double value, const CoeffsMatrix& coeffsmatrix) {
     445           0 :   return coeffsmatrix+value;
     446             : }
     447             : 
     448             : 
     449           0 : CoeffsMatrix operator+(const CoeffsMatrix& coeffsmatrix, const double value) {
     450           0 :   return CoeffsMatrix(coeffsmatrix)+=value;
     451             : }
     452             : 
     453             : 
     454           0 : CoeffsMatrix& CoeffsMatrix::operator+=(const std::vector<double>& values) {
     455           0 :   addToValues(values);
     456           0 :   return *this;
     457             : }
     458             : 
     459             : 
     460           0 : CoeffsMatrix operator+(const std::vector<double>& values, const CoeffsMatrix& coeffsmatrix) {
     461           0 :   return coeffsmatrix+values;
     462             : }
     463             : 
     464             : 
     465           0 : CoeffsMatrix operator+(const CoeffsMatrix& coeffsmatrix, const std::vector<double>& values) {
     466           0 :   return CoeffsMatrix(coeffsmatrix)+=values;
     467             : }
     468             : 
     469             : 
     470           0 : CoeffsMatrix& CoeffsMatrix::operator-=(const double value) {
     471           0 :   subtractFromValues(value);
     472           0 :   return *this;
     473             : }
     474             : 
     475             : 
     476           0 : CoeffsMatrix operator-(const double value, const CoeffsMatrix& coeffsmatrix) {
     477           0 :   return -1.0*coeffsmatrix+value;
     478             : }
     479             : 
     480             : 
     481           0 : CoeffsMatrix operator-(const CoeffsMatrix& coeffsmatrix, const double value) {
     482           0 :   return CoeffsMatrix(coeffsmatrix)-=value;
     483             : }
     484             : 
     485             : 
     486           0 : CoeffsMatrix& CoeffsMatrix::operator-=(const std::vector<double>& values) {
     487           0 :   subtractFromValues(values);
     488           0 :   return *this;
     489             : }
     490             : 
     491             : 
     492           0 : CoeffsMatrix operator-(const std::vector<double>& values, const CoeffsMatrix& coeffsmatrix) {
     493           0 :   return -1.0*coeffsmatrix+values;
     494             : }
     495             : 
     496             : 
     497           0 : CoeffsMatrix operator-(const CoeffsMatrix& coeffsmatrix, const std::vector<double>& values) {
     498           0 :   return CoeffsMatrix(coeffsmatrix)-=values;
     499             : }
     500             : 
     501             : 
     502           0 : CoeffsMatrix& CoeffsMatrix::operator+=(const CoeffsMatrix& other_coeffsmatrix) {
     503           0 :   addToValues(other_coeffsmatrix);
     504           0 :   return *this;
     505             : }
     506             : 
     507             : 
     508           0 : CoeffsMatrix CoeffsMatrix::operator+(const CoeffsMatrix& other_coeffsmatrix) const {
     509           0 :   return CoeffsMatrix(*this)+=other_coeffsmatrix;
     510             : }
     511             : 
     512             : 
     513           0 : CoeffsMatrix& CoeffsMatrix::operator-=(const CoeffsMatrix& other_coeffsmatrix) {
     514           0 :   subtractFromValues(other_coeffsmatrix);
     515           0 :   return *this;
     516             : }
     517             : 
     518             : 
     519           0 : CoeffsMatrix CoeffsMatrix::operator-(const CoeffsMatrix& other_coeffsmatrix) const {
     520           0 :   return CoeffsMatrix(*this)-=other_coeffsmatrix;
     521             : }
     522             : 
     523             : 
     524           0 : void CoeffsMatrix::averageMatrices(CoeffsMatrix& coeffsmat0, CoeffsMatrix& coeffsmat1) {
     525           0 :   plumed_massert(sameShape(coeffsmat0,coeffsmat1),"both CoeffsMatrix objects need to have the same shape");
     526           0 :   for(size_t i=0; i<coeffsmat0.getSize(); i++) {
     527           0 :     coeffsmat0.data[i] = coeffsmat1.data[i] = 0.5 * (coeffsmat0.data[i]+coeffsmat1.data[i]);
     528             :   }
     529           0 : }
     530             : 
     531             : 
     532           0 : void CoeffsMatrix::averageMatrices(const std::vector<CoeffsMatrix*>& coeffsmatSet) {
     533           0 :   double norm_factor = 1.0/static_cast<double>(coeffsmatSet.size());
     534           0 :   for(unsigned int k=1; k<coeffsmatSet.size(); k++) {
     535           0 :     plumed_massert(coeffsmatSet[0]->sameShape(*coeffsmatSet[k]),"All CoeffsMatrix objects need to have the same shape");
     536             :   }
     537           0 :   for(size_t i=0; i<coeffsmatSet[0]->getSize(); i++) {
     538             :     double value = 0.0;
     539           0 :     for(unsigned int k=0; k<coeffsmatSet.size(); k++) {
     540           0 :       value += coeffsmatSet[k]->data[i];
     541             :     }
     542           0 :     value *= norm_factor;
     543           0 :     for(unsigned int k=0; k<coeffsmatSet.size(); k++) {
     544           0 :       coeffsmatSet[k]->data[i] = value;
     545             :     }
     546             :   }
     547           0 : }
     548             : 
     549             : 
     550             : 
     551           0 : double CoeffsMatrix::getMinValue() const {
     552             :   double min_value=DBL_MAX;
     553           0 :   for(size_t i=0; i<data.size(); i++) {
     554           0 :     if(data[i]<min_value) {
     555             :       min_value=data[i];
     556             :     }
     557             :   }
     558           0 :   return min_value;
     559             : }
     560             : 
     561             : 
     562           0 : double CoeffsMatrix::getMaxValue() const {
     563             :   double max_value=DBL_MIN;
     564           0 :   for(size_t i=0; i<data.size(); i++) {
     565           0 :     if(data[i]>max_value) {
     566             :       max_value=data[i];
     567             :     }
     568             :   }
     569           0 :   return max_value;
     570             : }
     571             : 
     572             : 
     573           0 : void CoeffsMatrix::randomizeValuesGaussian(int randomSeed) {
     574           0 :   Random rnd;
     575             :   if (randomSeed<0) {randomSeed = -randomSeed;}
     576           0 :   rnd.setSeed(-randomSeed);
     577           0 :   for(size_t i=0; i<data.size(); i++) {
     578           0 :     data[i]=rnd.Gaussian();
     579             :   }
     580           0 : }
     581             : 
     582             : 
     583           0 : void CoeffsMatrix::resetAveraging() {
     584           0 :   clear();
     585             :   resetAveragingCounter();
     586           0 : }
     587             : 
     588             : 
     589           0 : void CoeffsMatrix::addToAverage(const CoeffsMatrix& coeffsmat) {
     590           0 :   plumed_massert( data.size()==coeffsmat.getSize(), "Incorrect size");
     591             :   //
     592           0 :   double aver_decay = 1.0 / ( static_cast<double>(averaging_counter) + 1.0 );
     593           0 :   if(averaging_exp_decay_>0 &&  (averaging_counter+1 > averaging_exp_decay_) ) {
     594           0 :     aver_decay = 1.0 / static_cast<double>(averaging_exp_decay_);
     595             :   }
     596             :   //
     597           0 :   for(size_t i=0; i<data.size(); i++) {
     598           0 :     data[i]+=(coeffsmat.data[i]-data[i])*aver_decay;
     599             :   }
     600           0 :   averaging_counter++;
     601           0 : }
     602             : 
     603             : 
     604         660 : void CoeffsMatrix::writeToFile(OFile& ofile) {
     605         660 :   writeHeaderToFile(ofile);
     606         660 :   writeDataToFile(ofile);
     607         660 : }
     608             : 
     609             : 
     610           0 : void CoeffsMatrix::writeToFile(const std::string& filepath, const bool append_file, Action* action_pntr) {
     611           0 :   OFile file;
     612           0 :   if(action_pntr!=NULL) {
     613           0 :     file.link(*action_pntr);
     614             :   }
     615             :   else {
     616           0 :     file.link(mycomm);
     617             :   }
     618           0 :   if(append_file) { file.enforceRestart(); }
     619           0 :   file.open(filepath);
     620           0 :   writeToFile(file);
     621           0 :   file.close();
     622           0 : }
     623             : 
     624             : 
     625         660 : void CoeffsMatrix::writeMatrixInfoToFile(OFile& ofile) {
     626         660 :   std::string field_diagonal = "diagonal_matrix";
     627         660 :   ofile.addConstantField(field_diagonal).printField(field_diagonal,isDiagonal());
     628         660 : }
     629             : 
     630             : 
     631         660 : void CoeffsMatrix::writeHeaderToFile(OFile& ofile) {
     632         660 :   ofile.clearFields();
     633         660 :   if(isIterationCounterActive()) {
     634         660 :     writeIterationCounterAndTimeToFile(ofile);
     635             :   }
     636         660 :   writeCoeffsInfoToFile(ofile);
     637         660 :   writeMatrixInfoToFile(ofile);
     638         660 : }
     639             : 
     640             : 
     641         660 : void CoeffsMatrix::writeDataToFile(OFile& ofile) {
     642         660 :   if(diagonal_) {
     643         660 :     writeDataDiagonalToFile(ofile);
     644             :   }
     645             :   else {
     646           0 :     writeDataFullToFile(ofile);
     647             :   }
     648         660 : }
     649             : 
     650             : 
     651         660 : void CoeffsMatrix::writeDataDiagonalToFile(OFile& ofile) {
     652             :   //
     653         660 :   std::string field_indices_prefix = "idx_";
     654             :   std::string field_coeffs = getDataLabel();
     655         660 :   std::string field_index = "index";
     656             :   //
     657         660 :   std::string int_fmt = "%8d";
     658         660 :   std::string str_separate = "#!-------------------";
     659             :   //
     660         660 :   std::vector<char> s1(20);
     661         660 :   std::vector<unsigned int> indices(numberOfDimensions());
     662         660 :   std::vector<std::string> ilabels(numberOfDimensions());
     663        1520 :   for(unsigned int k=0; k<numberOfDimensions(); k++) {
     664        1720 :     ilabels[k]=field_indices_prefix+getDimensionLabel(k);
     665             :   }
     666             :   //
     667       28290 :   for(size_t i=0; i<numberOfCoeffs(); i++) {
     668       27630 :     indices=getIndices(i);
     669       77540 :     for(unsigned int k=0; k<numberOfDimensions(); k++) {
     670       49910 :       std::snprintf(s1.data(),s1.size(),int_fmt.c_str(),indices[k]);
     671       99820 :       ofile.printField(ilabels[k],s1.data());
     672             :     }
     673       55260 :     ofile.fmtField(" "+getOutputFmt()).printField(field_coeffs,getValue(i,i));
     674       27630 :     std::snprintf(s1.data(),s1.size(),int_fmt.c_str(),i); ofile.printField(field_index,s1.data());
     675       27630 :     ofile.printField();
     676             :   }
     677         660 :   ofile.fmtField();
     678             :   // blank line between iterations to allow proper plotting with gnuplot
     679         660 :   ofile.printf("%s\n",str_separate.c_str());
     680         660 :   ofile.printf("\n");
     681         660 :   ofile.printf("\n");
     682        1320 : }
     683             : 
     684             : 
     685           0 : void CoeffsMatrix::writeDataFullToFile(OFile& ofile) {
     686             :   //
     687           0 :   std::string field_index_row = "idx_row";
     688           0 :   std::string field_index_column = "idx_column";
     689             :   std::string field_coeffs = getDataLabel();
     690             :   //
     691           0 :   std::string int_fmt = "%8d";
     692           0 :   std::string str_separate = "#!-------------------";
     693             :   //
     694           0 :   std::vector<char> s1(20);
     695             :   //
     696           0 :   for(size_t i=0; i<nrows_; i++) {
     697           0 :     for(size_t j=0; j<ncolumns_; j++) {
     698             :       std::snprintf(s1.data(),s1.size(),int_fmt.c_str(),i);
     699           0 :       ofile.printField(field_index_row,s1.data());
     700             :       std::snprintf(s1.data(),s1.size(),int_fmt.c_str(),j);
     701           0 :       ofile.printField(field_index_column,s1.data());
     702           0 :       ofile.fmtField(" "+getOutputFmt()).printField(field_coeffs,getValue(i,j));
     703           0 :       ofile.printField();
     704             :     }
     705             :   }
     706           0 :   ofile.fmtField();
     707             :   // blank line between iterations to allow proper plotting with gnuplot
     708           0 :   ofile.printf("%s\n",str_separate.c_str());
     709           0 :   ofile.printf("\n");
     710           0 :   ofile.printf("\n");
     711           0 : }
     712             : 
     713             : 
     714             : }
     715             : }