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