Kernel.h

Go to the documentation of this file.

/********************************************************************
 *  Copyright (C) 2010 by Federico Marulli                          *
 *  federico.marulli3@unibo.it                                      *
 *                                                                  *
 *  This program is free software; you can redistribute it and/or   *
 *  modify it under the terms of the GNU General Public License as  *
 *  published by the Free Software Foundation; either version 2 of  *
 *  the License, or (at your option) any later version.             *
 *                                                                  *
 *  This program is distributed in the hope that it will be useful, *
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of  *
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the   *
 *  GNU General Public License for more details.                    *
 *                                                                  *
 *  You should have received a copy of the GNU General Public       *
 *  License along with this program; if not, write to the Free      *
 *  Software Foundation, Inc.,                                      *
 *  59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.       *
 ********************************************************************/
 
#ifndef __KERNEL__
#define __KERNEL__
 
#include "Path.h"
#include "Constants.h"
#include "EnumCast.h"
#include "Exception.h"
 
// ============================================================================================
 
 
namespace cbl {
  
  enum class Dim {
    
          _1D_,
    
          _2D_
    
  };
 
  inline std::vector<std::string> DimNames () { return {"1D", "2D"}; }
  
  enum class BinType { 
 
              _linear_,
      
              _logarithmic_,
 
              _custom_
      
  };
 
  inline std::vector<std::string> BinTypeNames () { return {"linear", "logarithmic"}; }
 
  inline BinType BinTypeCast (const int binTypeIndex) { return castFromValue<BinType>(binTypeIndex); }
 
  inline BinType BinTypeCast (const std::string binTypeName) { return castFromName<BinType>(binTypeName, BinTypeNames()); }
 
  inline std::vector<BinType> BinTypeCast (const std::vector<int> binTypeIndeces) { return castFromValues<BinType>(binTypeIndeces); } 
 
  inline std::vector<BinType> BinTypeCast (const std::vector<std::string> binTypeNames) { return castFromNames<BinType>(binTypeNames, BinTypeNames()); }
 
  enum class CoordinateUnits {
 
                  _radians_,
    
                  _degrees_,
 
                  _arcseconds_,
 
                  _arcminutes_
    
  };
  
 
  inline std::vector<std::string> CoordinateUnitsNames () { return {"radians", "degrees", "arcseconds", "arcminutes"}; }
 
  inline CoordinateUnits CoordinateUnitsCast (const int coordinateUnitsIndex) { return castFromValue<CoordinateUnits>(coordinateUnitsIndex); }
 
  inline CoordinateUnits CoordinateUnitsCast (const std::string coordinateUnitsName) { return castFromName<CoordinateUnits>(coordinateUnitsName, CoordinateUnitsNames()); }
 
  inline std::vector<CoordinateUnits> CoordinateUnitsCast (const std::vector<int> coordinateUnitsIndeces) { return castFromValues<CoordinateUnits>(coordinateUnitsIndeces); } 
 
  inline std::vector<CoordinateUnits> CoordinateUnitsCast (const std::vector<std::string> coordinateUnitsNames) { return castFromNames<CoordinateUnits>(coordinateUnitsNames, CoordinateUnitsNames()); }
 
  
  enum class CoordinateType {
 
                 _comoving_,
    
                 _observed_
    
  };
 
  inline std::vector<std::string> CoordinateTypeNames () { return {"comoving", "observed"}; }
 
  inline CoordinateType CoordinateTypeCast (const int coordinateTypeIndex) { return castFromValue<CoordinateType>(coordinateTypeIndex); }
 
  inline CoordinateType CoordinateTypeCast (const std::string coordinateTypeName) { return castFromName<CoordinateType>(coordinateTypeName, CoordinateTypeNames()); }
 
  inline std::vector<CoordinateType> CoordinateTypeCast (const std::vector<int> coordinateTypeIndeces) { return castFromValues<CoordinateType>(coordinateTypeIndeces); } 
 
  inline std::vector<CoordinateType> CoordinateTypeCast (const std::vector<std::string> coordinateTypeNames) { return castFromNames<CoordinateType>(coordinateTypeNames, CoordinateTypeNames()); }
 
  struct comovingCoordinates { double xx; double yy; double zz; };
  struct observedCoordinates { double ra; double dec; double redshift; };
 
  typedef Eigen::Matrix<double, 3, 1> Vector3D;
 
  typedef Eigen::Matrix<double, 4, 1> Vector4D;
 
  typedef Eigen::Matrix<std::complex<double>, 1, Eigen::Dynamic> VectorComplex;
 
  typedef std::function<double(double)> FunctionDoubleDouble;
 
  typedef std::function<double(std::vector<double>)> FunctionDoubleVector;
 
  typedef std::function<double(std::vector<double> &)> FunctionDoubleVectorRef;
 
  typedef std::function< double(double, std::shared_ptr<void>, std::vector<double> &)> FunctionDoubleDoublePtrVectorRef;
 
  typedef std::function< double(double, double, std::shared_ptr<void>, std::vector<double> &)> FunctionDoubleDoubleDoublePtrVectorRef;
 
  typedef std::function< double(std::vector<double>, std::shared_ptr<void>, std::vector<double> &)> FunctionDoubleVectorPtrVectorRef;
 
  typedef std::function< std::vector<double>(std::vector<double>, std::shared_ptr<void>, std::vector<double> &)> FunctionVectorVectorPtrVectorRef;
 
  typedef std::vector<std::vector<std::vector<int>>> Tensor3Di;
    
  typedef std::vector<std::vector<std::vector<double>>> Tensor3Dd;
    
  typedef std::vector<std::vector<std::vector<std::vector<int>>>> Tensor4Di;
         
 
 
  inline std::ostream &headerCBL (std::ostream &stream)
  {
    stream << par::col_blue << "CBL > " << par::col_default;
    return stream;
  }
 
#define coutCBL std::cout << headerCBL
  
  inline void WarningMsgCBL (const std::string msg, const std::string functionCBL, const std::string fileCBL)
  { std::cerr << std::endl << par::col_bred << "CBL > Warning in the CBL function " << cbl::par::col_yellow << functionCBL << cbl::par::col_bred << " of " << fileCBL << ": " << cbl::par::col_default << msg << std::endl << std::endl; }
 
  inline int Error (const std::string msg, const cbl::glob::ExitCode exitCode=cbl::glob::ExitCode::_error_, const std::string header="\n")
  { throw cbl::glob::Exception(msg, exitCode, header); }
 
  inline int ErrorCBL (const std::string msg, const std::string functionCBL, const std::string fileCBL, const cbl::glob::ExitCode exitCode=cbl::glob::ExitCode::_error_)
  { throw cbl::glob::Exception(msg, exitCode, cbl::par::ErrorMsg, functionCBL, fileCBL); }
  
  inline void Beep (const std::string beep="beep")
  { if (system(("command -v say >/dev/null && { say "+beep+"; }").c_str())) {} }
 
  inline bool isSet (const std::string var) 
  { return (var!=cbl::par::defaultString) ? true : false; }
 
  inline bool isSet (const int var) 
  { return (var>par::defaultInt) ? true : false; }
 
  inline bool isSet (const long var) 
  { return (var>par::defaultLong) ? true : false; }
  
  inline bool isSet (const double var) 
  { return (var>par::defaultDouble) ? true : false; }
  
  inline bool isSet (const std::vector<double> vect) 
  {
    bool is = true;
    size_t ind = 0;
    while (is && ind<vect.size()) 
      if (vect[ind++]<par::defaultDouble*1.000001) is = false;
    return is;
  }
  
  inline bool isSet (const std::vector<unsigned int> vect) 
  {
    bool is = true;
    size_t ind = 0;
    while (is && ind<vect.size()) 
      if (vect[ind++]<par::defaultInt*1.000001) is = false;
    return is;
  }
  
  inline bool isSet (const std::vector<std::vector<unsigned int>> vect) 
  {
    bool is = true;
    size_t ind = 0;
    size_t i=0;
    while (is && i<vect.size()) {
        while (is && ind<vect[i].size()) 
            if (vect[i][ind++]<par::defaultInt*1.000001) is = false;
        i++;
    }
    return is;
  }
 
  template <typename T> std::string conv (const T val, const char *fact)
  {
    char VAL[20]; sprintf(VAL, fact, val); 
    return std::string(VAL);
  }
  
  template <typename T> 
  int nint (const T val) 
  { return (val<0) ? val-0.5 : val+0.5; }
  
  template <typename T> 
  T Log (const T val, const double fact=0.9) 
  { return (val>0) ? log10(val) : par::defaultDouble*fact; }
  
  template <typename T> 
  T Ln (const T val, const double fact=0.9) 
  { return (val>0) ? log(val) : par::defaultDouble*fact; }
 
  template <typename T>
  T closest (T x, T a, T b)
  { 
    if (a>b) ErrorCBL("the input parameter a must be <= than the input parameter b!", "closest", "Kernel.h");
    else if (a==b) return a;
    else return (fabs(x-a) < fabs(x-b)) ? a : b;
    return 1;
  }
 
  template <typename T>
  T index_closest (T x, std::vector<T> vv)
  { 
    if (vv.size()==0) ErrorCBL("vv is an empty std::vector!", "index_closest", "Kernel.h");
    std::vector<double>::iterator low, up;
    low = lower_bound(vv.begin(), vv.end(), x);
    up = upper_bound(vv.begin(), vv.end(), x);
    int index = (closest(x, *low, *up)==*low) ? low-vv.begin() : up-vv.begin();
    return index;
  }
 
  template <typename T>
  T closest (T x, std::vector<T> values)
  { return values[index_closest(x, values)]; }
  
  short ShortSwap (const short s);
 
  int IntSwap (const int i);
 
  long long LongSwap (const long long i);
 
  float FloatSwap (const float f);
 
  double DoubleSwap (const double d);
 
  double round_to_digits (const double num, const int ndigits);
 
  double round_to_precision (const double num, const int ndigits);
 
  void checkIO (const std::ifstream &fin, const std::string file="NULL");
 
  void checkIO (const std::ofstream &fout, const std::string file="NULL");
 
  void set_EnvVar (const std::vector<std::string> Var);
 
  void check_EnvVar (const std::string Var); 
 
  int used_memory (const int type);
 
  int check_memory (const double frac, const bool exit=true, const std::string func="", const int type=1);
 
 
  
  // ============================================================================================
  
 
 
  template <typename T> 
  void Print (const T value, const int prec, const int ww, const std::string header="", const std::string end="\n", const bool use_coutCBL=true, std::ostream &stream=std::cout, const std::string colour=cbl::par::col_default) 
  {
    const int bp = std::cout.precision(); 
    if (fabs(value)<pow(10, -prec) || fabs(value)>pow(10, prec)) {
      if (use_coutCBL)
    coutCBL << header << colour << std::scientific << std::setprecision(prec) << std::setw(ww) << std::right << value << par::col_default << end;
      else
    stream << header << std::scientific << std::setprecision(prec) << std::setw(ww) << std::right << value << end;
    }
    else {
      if (use_coutCBL) 
    coutCBL << header << colour << std::fixed << std::setprecision(prec) << std::setw(ww) << std::right << value << par::col_default << end;
      else 
    stream << header << std::fixed << std::setprecision(prec) << std::setw(ww) << std::right << value << end;
    }
    std::cout.precision(bp);
  }
  
  template <typename T> 
  void Print (const std::vector<T> vect, const int prec=4, const int ww=8) 
  {
    for (auto &&vi : vect)
      Print(vi, prec, ww);
  }
  
  inline void Print (const std::vector<std::string> vect)
  {
    for (auto &&vi : vect) 
      coutCBL << vi << std::endl;
  }
  
 
  template <typename T> 
  void Print (const std::vector<T> vect1, const std::vector<T> vect2, const int prec=4, const int ww=8) 
  {
    if (vect1.size()!=vect2.size())
      ErrorCBL("the two input vectors to be printed must have the same dimension!", "Print", "Kernel.h");
      
    for (size_t i=0; i<vect1.size(); i++) {
      Print(vect1[i], prec, ww, "", "  ");
      Print(vect2[i], prec, ww);
    }
  }
 
 
  inline void Print (const std::vector<std::string> vect1, const std::vector<std::string> vect2) 
  {
    if (vect1.size()!=vect2.size())
      ErrorCBL("the two input vectors to be printed must have the same dimension!", "Print", "Kernel.h");
      
    for (size_t i=0; i<vect1.size(); i++) 
      coutCBL << vect1[i] << "   " << vect2[i] << std::endl;
  }
 
 
  template <typename T> 
  void Print (const std::vector<T> vect1, const std::vector<T> vect2, const std::vector<T> vect3, const int prec=4, const int ww=8) 
  {
    if (vect1.size()!=vect2.size() || vect1.size()!=vect3.size())
      ErrorCBL("the three input vectors to be printed must have the same dimension!", "Print", "Kernel.h");
      
    for (size_t i=0; i<vect1.size(); i++) {
      Print(vect1[i], prec, ww, "", "  ");
      Print(vect2[i], prec, ww, "", "  ");
      Print(vect3[i], prec, ww);    
    }
  }
 
 
  inline void Print (const std::vector<std::string> vect1, const std::vector<std::string> vect2, const std::vector<std::string> vect3) 
  {
    if (vect1.size()!=vect2.size() || vect1.size()!=vect3.size())
      ErrorCBL("the three input vectors to be printed must have the same dimension!", "Print", "Kernel.h");
      
    for (size_t i=0; i<vect1.size(); i++) 
      coutCBL << vect1[i] << "   " << vect2[i] << "   " << vect3[i] << std::endl;
  }
  
  
  template <typename T> 
  void Print (const std::vector<std::vector<T>> mat, const int prec=4, const int ww=8) 
  {
    for (size_t i=0; i<mat.size(); i++) {
      for (size_t j=0; j<mat[i].size(); j++)
    if (j==0)
      Print(mat[i][j], prec, ww, "", "  ");
    else
      Print(mat[i][j], prec, ww, "", "  "); 
      std::cout << std::endl;
    }
  }
 
  inline void Print (const std::vector<std::vector<std::string>> mat) 
  {
    for (size_t i=0; i<mat.size(); i++) {
      for (size_t j=0; j<mat[i].size(); j++)
    if (j==0) 
      coutCBL << mat[i][j] << "   ";
    else 
      std::cout << mat[i][j] << "   ";
      std::cout << std::endl;
    }
  }
  
 
  template <typename T> 
  T Min (const std::vector<T> vect) 
  {
    if (vect.size()==0) ErrorCBL("vect.size=0!", "Min", "Kernel.h");
    return *min_element(vect.begin(), vect.end());
  }
 
  template <typename T> 
  T Max (const std::vector<T> vect) 
  {
    if (vect.size()==0) ErrorCBL("vect.size=0!", "Max", "Kernel.h");
    return *max_element(vect.begin(), vect.end());
  }
 
  template <typename T> 
  std::vector<T> different_elements (const std::vector<T> vect_input) 
  {
    std::vector<T> vect = vect_input;
    sort(vect.begin(), vect.end());
    typename std::vector<T>::iterator it = unique(vect.begin(), vect.end()); 
    vect.resize(it-vect.begin());    
    return vect;
  }
 
  template <typename T> 
  int N_different_elements (const std::vector<T> vect_input) 
  {
    std::vector<T> vect = different_elements<T>(vect_input);
    return vect.size();
  }
 
  void unique_unsorted (std::vector<int> & vv);
 
  void unique_unsorted (std::vector<double> & vv);
 
  void unique_unsorted (std::vector<std::string> & vv);
 
  template <typename T> 
  void Erase (std::vector<T> &vv, std::vector<int> ind) 
  {
    for (auto &&i : ind) 
      if (i>=int(vv.size())) ErrorCBL("the input value of ind is too large!", "Erase", "Kernel.h");
 
    unique_unsorted(ind);
    int tt = 0;
    for (auto &&i : ind) 
      vv.erase(vv.begin()+i-(tt++));
  }
 
  template <typename T> 
  void Erase_lines (std::vector<std::vector<T> > &Mat, std::vector<int> ll) 
  {
    for (auto &&i : ll)
      if (i>=int(Mat.size())) ErrorCBL("the dimension of the input vector ll is too large!", "Erase_lines", "Kernel.h");
 
    unique_unsorted(ll);
    int tt = 0;
    for (auto &&i : ll)
      Mat.erase(Mat.begin()+i-(tt++));
  }
 
  template <typename T> 
  void Erase_columns (std::vector<std::vector<T> > &Mat, std::vector<int> col) 
  {
    for (auto &&i : col)
      for (auto &&j : Mat)
    if (i>=int(j.size())) ErrorCBL("the dimension of the input vector col is too large!", "Erase_columns", "Kernel.h");
 
    unique_unsorted(col);
    int tt = 0;
    for (auto &&i : col) {
      for (auto &&j : Mat)
    j.erase(j.begin()+i-tt);
      tt ++;
    }
  }
 
  template <typename T> 
  void SubMatrix (std::vector<T> &xx, std::vector<T> &yy, std::vector<std::vector<T> > &Mat, T val) 
  { 
    std::vector<int> line, column;
 
    for (unsigned int i=0; i<xx.size(); i++) {
      if (i>=Mat.size()) ErrorCBL("the dimension of the input vector xx is too large!", "SubMatrix", "Kernel.h");
      bool ll = 0;
 
      for (unsigned int j=0; j<yy.size(); j++) {
    if (j>=Mat[i].size()) ErrorCBL("the dimension of the input vector yy is too large!", "SubMatrix", "Kernel.h");
    if (Mat[i][j]<val) {
      if (j<int(yy.size()*0.5)) {line.push_back(i); ll = 1;}
      else if (ll==0) {column.push_back(j);}
    }
      }
    }
      
    Erase(xx, line);
    Erase_lines(Mat, line);
      
    Erase(yy, column);
    Erase_columns(Mat, column);
  }
 
  template <typename T> 
  bool isDimEqual (const std::vector<T> vect1, const std::vector<T> vect2) 
  {
    return (vect1.size()==vect2.size()) ? 1 : 0;
  }
 
  template <typename T> 
  bool isDimEqual (const std::vector<std::vector<T> > mat1, const std::vector<std::vector<T> > mat2) 
  {
    bool is = (mat1.size()==mat2.size()) ? 1 : 0;
    if (is) 
      for (unsigned int i=0; i<mat1.size(); i++) {
    if (mat1[i].size()!=mat2[i].size()) is = 0;
      }
    return is;
  }
 
  template <typename T> 
  void checkDim (const std::vector<T> vect, const int val, const std::string vector, bool equal=true) 
  {
    if (equal) {
      if ((int)vect.size()!=val) 
    ErrorCBL("the dimension of " + vector + " is: " + conv(vect.size(), par::fINT) + " ( != " + conv(val, par::fINT) + " )", "checkDim", "Kernel.h");
    }
    else { 
      if ((int)vect.size()<val)
    ErrorCBL("the dimension of " + vector + " is: " + conv(vect.size(), par::fINT) + " ( < " + conv(val, par::fINT) + " )", "checkDim", "Kernel.h");
    }
  }
  
  template <typename T> 
  void checkDim (const std::vector<T> mat, const int val_i, const int val_j, const std::string matrix, const bool equal=true) 
  {
    if (equal) {
      if (int(mat.size())!=val_i) 
    ErrorCBL("the dimension of " + matrix + " is: " + conv(mat.size(), par::fINT) + " != " + conv(val_i, par::fINT) + "!", "checkDim", "Kernel.h");
      else 
    for (size_t k=0; k<mat.size(); k++)
      if (int(mat[k].size())!=val_j) 
        ErrorCBL("the dimension of " + matrix + " is: " + conv(mat[k].size(), par::fINT) + " != " + conv(val_j, par::fINT) + "!", "checkDim", "Kernel.h");
    }
    else {
      if (int(mat.size())<val_i) 
    ErrorCBL("the dimension of " + matrix + " is: " + conv(mat.size(), par::fINT) + " < " + conv(val_i, par::fINT) + "!", "checkDim", "Kernel.h");
      else 
    for (size_t k=0; k<mat.size(); k++)
      if (int(mat[k].size())<val_j) 
        ErrorCBL("the dimension of " + matrix + " is: " + conv(mat[k].size(), par::fINT) + " < " + conv(val_j, par::fINT) + "!", "checkDim", "Kernel.h");
    }
  }
 
  template <typename T> 
  void checkEqual (const std::vector<T> vect1, const std::vector<T> vect2) 
  {
    checkDim(vect2, vect1.size(), "vect2");
    for (size_t i=0; i<vect1.size(); i++)
      if (vect1[i]!=vect2[i])
    ErrorCBL("vect1 and vect2 are different!", "checkEqual", "Kernel.h");
  }
  
  template <typename T> 
  std::vector<T> linear_bin_vector (const size_t nn, const T min, const T max)
  {
    std::vector<T> vv(nn);
    for (size_t i = 0; i<nn; i++)
      vv[i] = min+(max-min)*T(i)/T(nn-1);
    return vv;
  }
 
  template <typename T> 
  std::vector<T> logarithmic_bin_vector (const size_t nn, const T min, const T max)
  {
    std::vector<T> vv(nn);
    for (size_t i=0; i<nn; i++)
      vv[i] = exp(log(min)+(log(max)-log(min))*T(i)/T(nn-1));
    return vv;
  }
 
  template <typename T> 
  int locate (const std::vector<T> &vv, const T xx) 
  {
    size_t nn = vv.size ();
    int jl = -1;
    int ju = nn;
    bool as = (vv[nn-1] >= vv[0]);
    while (ju-jl > 1)
      {
    int jm = (ju+jl)*0.5;
    if ((xx >= vv[jm]) == as)
      jl = jm;
    else
      ju = jm;
      }
    if (xx == vv[0])
      return 0;
    else if (xx == vv[nn-1])
      return nn-2;
    else
      return jl;
  }
 
  template <typename T>
  std::vector<T> extract_elements (std::vector<T> vec, std::vector<unsigned int> index)
  {
    std::vector<T> vv;
    for (unsigned int i=0; i< index.size(); i++)
      vv.push_back(vec[index[i]]);
    return vv;
  }
  
  template <typename T>
  std::vector<T> flatten(std::vector<std::vector<T>> matrix)
  {
    std::vector<T> flatten;
    for (size_t i=0; i<matrix.size(); i++)
      for (size_t j=0; j<matrix[i].size(); j++)
    flatten.push_back(matrix[i][j]);
 
    return flatten;
  }
 
  template <typename T>
  std::vector<std::vector<T>> reshape (std::vector<T> vec, const int size1, const int size2)
  {
    if (size1*size2!=int(vec.size()))
      ErrorCBL("sizes does not match! "+conv(size1*size2, par::fINT)+" should be equal to "+conv(int(vec.size()), par::fINT), "reshape", "Kernel.h");
 
    std::vector<std::vector<T>> matrix(size1, std::vector<T> (size2, 0));
 
    for (int i=0; i<size1; i++)
      for (int j=0; j<size2; j++)
    matrix[i][j] = vec[j+i*size2];
 
    return matrix;
  }
  
  template <typename T>
  std::vector<std::vector<T>> transpose (std::vector<std::vector<T>> matrix)
  {
    const int size1 = matrix.size();
    const int size2 = matrix[0].size();
      
    std::vector<std::vector<T>> TRmatrix(size2, std::vector<T> (size1, 0));
 
    for (int i=0; i<size1; i++)
      for (int j=0; j<size2; j++)
    TRmatrix[j][i] = matrix[i][j];
 
    return TRmatrix;
  }
 
  template <typename T>
  bool inRange (T value, T min, T max, bool include_limits = true) 
  {
    if (include_limits){
      if (value>=min && max>=value)
    return true;
    }
    else{
      if (value>min && max>value)
    return true;
    }
 
    return false;
  }
 
  template <typename T>
  bool inRange (std::vector<T> value, std::vector<T> min, std::vector<T> max, bool include_limits = true) 
  {
    int in_range = 1;
 
    for (size_t i=0; i<value.size(); i++)
      in_range *= int(inRange(value[i], min[i], max[i], include_limits));
 
    return bool(in_range);
  }
 
  template <typename T>
  bool inRange (std::vector<T> value, std::vector<std::vector<T>> ranges, bool include_limits = true) 
  {
    int in_range = 1;
 
    for (size_t i=0; i<value.size(); i++)
      in_range *= int(inRange(value[i], ranges[i][0], ranges[i][1], include_limits));
 
    return bool(in_range);
  }
 
  template <typename T> 
  T v_M_vt (const std::vector<T> vv, const std::vector<std::vector<T>> MM)
  {
    const int size = vv.size();
 
    std::vector<double> ivv(size, 0);
    for(int i=0; i<size; i++)
      for(int j=0; j<size; j++)
    ivv[i] += vv[j]*MM[i][j];
 
    double res=0;
    for(int i=0; i<size; i++)
      res += vv[i]*ivv[i];
 
    return res;
  }
 
  // sort two or more std::vectors at the same time
  namespace glob {
    class CL {
    public:
      std::vector<double> VV;
      CL (std::vector<double> vv) {VV = vv;};
    };
    bool operator<(const CL &, const CL &);
  }
  
  void sort_2vectors (std::vector<double>::iterator p1, std::vector<double>::iterator p2, const int dim);
 
  void sort_3vectors (std::vector<double>::iterator p1, std::vector<double>::iterator p2, std::vector<double>::iterator p3, const int dim);
 
  void sort_4vectors (std::vector<double>::iterator p1, std::vector<double>::iterator p2, std::vector<double>::iterator p3, std::vector<double>::iterator p4, const int dim);
 
  int makeDir (std::string path, const std::string rootPath=".", const mode_t mode=0777, const bool verbose=false);
  
  inline std::vector<std::vector<double> > operator * (const std::vector<std::vector<double> > &Mat1, const std::vector<std::vector<double> > &Mat2)
  {   
    std::vector<std::vector<double> > MatP(Mat1.size(), std::vector<double>(Mat2[0].size(),0.));
  
    for (unsigned int i=0; i<Mat1.size(); i++) 
      for (unsigned int j=0; j<Mat2[0].size(); j++) {
    double temp = 0.;
    for (unsigned int k=0; k<Mat1[0].size(); k++) 
      temp += Mat1[i][k]*Mat2[k][j];
    MatP[i][j] = temp;
      }
 
    return MatP;
  }
 
  template <typename T>
  std::vector<T> slice (const std::vector<T> v, const int start=0, const  int end=-1)
  {
    int oldlen = v.size();
    int newlen;
 
    if (end==-1 || end>=oldlen)
      newlen = oldlen-start;
    else 
      newlen = end-start;
      
    std::vector<T> nv(newlen);
 
    for (int i=0; i<newlen; i++) 
      nv[i] = v[start+i];
      
    return nv;
  }
  
  std::vector<std::vector<double>> read_file (const std::string file_name, const std::string path_name, const std::vector<int> column_data, const int skip_nlines=0);
  
  std::vector<std::vector<double>> read_file (const std::string file_name, const std::string path_name, const std::vector<int> column_data, const std::string delimiter, const char comment='#');
  
}
 
 
#endif