MC_double_vector.cpp

Go to the documentation of this file.

#include <cmath>
#include <MC_double_vector.hpp>
#include <MC_string_helper.hpp>

#include <MC_int_vector.hpp>
#include <MC_v3d.hpp>
#include <MC_v3d_vector.hpp>

//using namespace mc;
namespace mesh_conv
{


  MC_double_vector::MC_double_vector(){}

  MC_double_vector::MC_double_vector(const double& u0)
  {double_list.push_back(u0);}
  MC_double_vector::MC_double_vector(const double& u0,const double& u1){double_list.push_back(u0);double_list.push_back(u1);}
  MC_double_vector::MC_double_vector(const double& u0,const double& u1,const double& u2){double_list.push_back(u0);double_list.push_back(u1);double_list.push_back(u2);}
  MC_double_vector::MC_double_vector(const double& u0,const double& u1,const double& u2,const double& u3){double_list.push_back(u0);double_list.push_back(u1);double_list.push_back(u2);double_list.push_back(u3);}
  MC_double_vector::MC_double_vector(const double& u0,const double& u1,const double& u2,const double& u3,const double& u4){double_list.push_back(u0);double_list.push_back(u1);double_list.push_back(u2);double_list.push_back(u3);double_list.push_back(u4);}
  MC_double_vector::MC_double_vector(const double& u0,const double& u1,const double& u2,const double& u3,const double& u4,const double& u5){double_list.push_back(u0);double_list.push_back(u1);double_list.push_back(u2);double_list.push_back(u3);double_list.push_back(u4);double_list.push_back(u5);}

  void MC_double_vector::assert_bounds(const int& u) const
  {
    if(u<0 || u>=int(double_list.size()))
      {std::cout<<"Assert MC_double_vector::assert_bounds("<<u<<") failed, because size="<<double_list.size()<<std::endl;exit(-1);}
  }  

  int MC_double_vector::size() const
  {return double_list.size();}

  const double& MC_double_vector::operator()(const int& index) const
  {assert_bounds(index); return double_list[index];}
  double& MC_double_vector::operator()(const int& index)
  {assert_bounds(index); return double_list[index];}
  const double& MC_double_vector::operator[](const int& index) const
  {assert_bounds(index); return double_list[index];}
  double& MC_double_vector::operator[](const int& index)
  {assert_bounds(index); return double_list[index];}

  void MC_double_vector::clear()
  {double_list.clear();}



  std::ostream& operator<<(std::ostream& output,const MC_double_vector& in)
  {
    int N=in.size();
    for(int k=0;k<N;k++)
      output<<in[k]<<" ";
    return output;
  }

  MC_double_vector::MC_double_vector(const std::string& in,const std::string& delimiter)
  {
    std::vector <std::string> v_string=MC_string_tokenizer::tokenize(in,delimiter);
    bool is_converted=false;
    for(unsigned int k=0;k<v_string.size();k++)
      {
        double value=MC_string_converter::value_of<double>(v_string[k],&is_converted);
        if(is_converted==true)
          add(value);
      }
  }  
  

  MC_double_vector& MC_double_vector::add(const double& val)
  {double_list.push_back(val);return *this;}

  MC_double_vector operator<<(const MC_double_vector& vec,const double& value)
  {return MC_double_vector(vec).add(value);}



  MC_double_vector& MC_double_vector::add(const MC_double_vector& val)
  {
      int N=val.size();
      int end=double_list.size();
      double_list.resize(double_list.size()+N);
      for(int k=0;k<N;++k)
          double_list[end+k]=val[k];
      return *this;
  }
  
  std::vector <double>& MC_double_vector::get_vector() {return double_list;}
  const std::vector <double>& MC_double_vector::to_vector() const {return double_list;}

  
  std::set <double> MC_double_vector::to_set() const
  {
    std::set <double> set_double;
    set_double.insert(double_list.begin(),double_list.end());
    return set_double;
  }
  std::map <double,int> MC_double_vector::to_map() const
  {
    std::map <double,int> map_double;
    int N=size();
    for(int k=0;k<N;k++)
      map_double.insert(std::pair<double,int>(this->operator()(k),k));
    return map_double;
  }

  MC_double_vector operator<<(const MC_double_vector& vec,const MC_double_vector& to_add)
  {return MC_double_vector(vec).add(to_add);}

  MC_double_vector& MC_double_vector::set(const int& k_index,const double& new_value)
  {
    if(k_index<0)
      {std::cout<<"Error in MC_double_vector::set("<<k_index<<","<<new_value<<"), index must be >0"<<std::endl;exit(-1);}
    if(k_index>=size())
      double_list.resize(k_index+1);
    double_list[k_index]=new_value;
    return *this;
  }
  
  MC_double_vector& MC_double_vector::set(const MC_int_vector& k_index,const MC_double_vector& new_value)
  {
    MC_double_vector v_new_value=new_value;
    if(new_value.size()!=1 && new_value.size()!=k_index.size())
      {std::cout<<"Error in MC_double_vector::set(MC_int_vector,MC_double_vector), size are not correct: "<<k_index.size()<<","<<new_value.size()<<std::endl;exit(-1);}
    if(new_value.size()==1)
      v_new_value=MC_double_vector::zeros(k_index.size())+new_value;
    
    
    int N=v_new_value.size();
    for(int k=0;k<N;k++)
      set(k_index[k],v_new_value[k]);
    return *this;
  }
    

  MC_double_vector::MC_double_vector(const std::set <double>& vec)
  {
    for(std::set <double> :: const_iterator it=vec.begin();it!=vec.end();++it)
        double_list.push_back(*it);
  }
  MC_double_vector::MC_double_vector(const std::map <double,int> &vec)
  {
    for(std::map <double,int> :: const_iterator it=vec.begin();it!=vec.end();++it)
      set(it->second,it->first);
  }
  MC_double_vector::MC_double_vector(const std::vector <double> &vec)
  {double_list=vec;}

  MC_double_vector::MC_double_vector(const MC_double_vector& vec)
  {double_list=vec.to_vector();}

  MC_double_vector MC_double_vector::zeros(const int& new_size)
  {
    MC_double_vector vec;vec.resize(new_size);
    return vec;
  }

  MC_double_vector& MC_double_vector::resize(const int& new_size)
  {
    int N=double_list.size();
    if(N==new_size) //same size
      return *this;
    
    if(new_size<0)
      {std::cout<<"Error in MC_double_vector::resize("<<new_size<<") new size must be >0"<<std::endl;exit(-1);}
    
    double_list.resize(new_size);
    return *this;
    
  }
  MC_double_vector operator-(const MC_double_vector& v)
  {
    MC_double_vector opposite;
    int N=v.size();
    opposite.resize(N);
    for(int k=0;k<N;++k)
      opposite[k]=-v[k];
    return opposite;
  }

  MC_double_vector operator+(const MC_double_vector& vec,const double& to_add)
  {
    int N=vec.size();
    MC_double_vector new_vec=MC_double_vector::zeros(N);
    for(int k=0;k<N;k++)
      new_vec[k]=vec[k]+to_add;
    return new_vec;
  }
  MC_double_vector operator-(const MC_double_vector& vec,const double& to_sub)
  {
    int N=vec.size();
    MC_double_vector new_vec=MC_double_vector::zeros(N);
    for(int k=0;k<N;k++)
      new_vec[k]=vec[k]-to_sub;
    return new_vec;
  }



  MC_double_vector operator*(const MC_double_vector& vec,const double& to_mult)
  {
      int N=vec.size();
      MC_double_vector new_vec=MC_double_vector::zeros(N);
      for(int k=0;k<N;k++)
          new_vec[k]=vec[k]*to_mult;
      return new_vec;
  }
  MC_double_vector operator/(const MC_double_vector& vec,const double& to_subdiv)
  {
      double epsilon=0.0000000001;
      if(fabs(to_subdiv)<epsilon)
      {std::cout<<"Error in MC_double_vector::operator/(MC_double_vector,double), double value is 0"<<std::endl;exit(-1);}
      int N=vec.size();
      MC_double_vector new_vec=MC_double_vector::zeros(N);
      for(int k=0;k<N;k++)
          new_vec[k]=vec[k]/to_subdiv;
      return new_vec;
  }

  MC_double_vector& MC_double_vector::operator+=(const double& to_add)
  {
    int N=size();
    for(int k=0;k<N;++k)
      double_list[k]+=to_add;
    return *this;
  }
  MC_double_vector& MC_double_vector::operator-=(const double& to_sub)
  {
    int N=size();
    for(int k=0;k<N;++k)
        double_list[k]-=to_sub;
    return *this;
  }
  MC_double_vector& MC_double_vector::operator*=(const double& to_mult)
  {
    int N=size();
    for(int k=0;k<N;++k)
      double_list[k]*=to_mult;
    return *this;
  }
  MC_double_vector& MC_double_vector::operator/=(const double& to_subdiv)
  {
      double epsilon=0.000001;
      if(fabs(to_subdiv)<epsilon)
      {std::cout<<"Error in MC_double_vector::operator/=(double), double value is 0"<<std::endl;exit(-1);}
      int N=size();
      for(int k=0;k<N;++k)
          double_list[k]/=to_subdiv;
      return *this;
  }


  MC_double_vector operator+(const MC_double_vector& vec,const MC_double_vector& to_add)
  {
    if(to_add.size()==1)
      return vec+to_add.first();
    else
      {
        if(vec.size()!=to_add.size())
          {std::cout<<"Error in MC_double_vector::operator+(MC_double_vector,MC_double_vector), size are not compatible "<<vec.size()<<";"<<to_add.size()<<std::endl;exit(-1);}
        

        int N=vec.size();
        MC_double_vector new_vec=MC_double_vector::zeros(N);
        for(int k=0;k<N;++k)
          new_vec[k]=vec[k]+to_add[k];
        return new_vec;
      }
  }
  MC_double_vector operator-(const MC_double_vector& vec,const MC_double_vector& to_sub)
  {
      if(to_sub.size()==1)
          return vec-to_sub.first();
      else
      {
          if(vec.size()!=to_sub.size())
          {std::cout<<"Error in MC_double_vector::operator-(MC_double_vector,MC_double_vector), size are not compatible "<<vec.size()<<";"<<to_sub.size()<<std::endl;exit(-1);}
        

          int N=vec.size();
          MC_double_vector new_vec=MC_double_vector::zeros(N);
          for(int k=0;k<N;++k)
              new_vec[k]=vec[k]-to_sub[k];
          return new_vec;
      }
  }
  MC_double_vector operator*(const MC_double_vector& vec,const MC_double_vector& to_mult)
  {
      if(to_mult.size()==1)
          return vec*to_mult.first();
      else
      {
          if(vec.size()!=to_mult.size())
          {std::cout<<"Error in MC_double_vector::operator*(MC_double_vector,MC_double_vector), size are not compatible "<<vec.size()<<";"<<to_mult.size()<<std::endl;exit(-1);}


          int N=vec.size();
          MC_double_vector new_vec=MC_double_vector::zeros(N);
          for(int k=0;k<N;++k)
              new_vec[k]=vec[k]*to_mult[k];
          return new_vec;
      }
  }
  MC_double_vector operator/(const MC_double_vector& vec,const MC_double_vector& to_subdiv)
  {
    if(to_subdiv.size()==1)
      return vec/to_subdiv.first();
    else
      {
        if(vec.size()!=to_subdiv.size())
          {std::cout<<"Error in MC_double_vector::operator/(MC_double_vector,MC_double_vector), size are not compatible "<<vec.size()<<";"<<to_subdiv.size()<<std::endl;exit(-1);}
        

        int N=vec.size();
        MC_double_vector new_vec=MC_double_vector::zeros(N);
        for(int k=0;k<N;++k)
          {
            double current=to_subdiv[k];
            double epsilon=0.0000001;
            if(fabs(current)<epsilon)
              {std::cout<<"Error in MC_double_vector::operator/(MC_double_vector,MC_double_vector), int value at index"<<k<<" is 0"<<std::endl;exit(-1);}
            new_vec[k]=vec[k]/current;
          }
        return new_vec;
      }

  }


  MC_double_vector& MC_double_vector::operator+=(const MC_double_vector& to_add)
  {
    if(to_add.size()==1)
      return *this+=to_add.first();
    else
      {
        if(size()!=to_add.size())
          {std::cout<<"Error in MC_double_vector::operator+=(MC_double_vector), size are not compatible "<<size()<<";"<<to_add.size()<<std::endl;exit(-1);}
        
        int N=size();
        for(int k=0;k<N;++k)
          double_list[k]+=to_add[k];
        return *this;
      }
  }
  MC_double_vector& MC_double_vector::operator-=(const MC_double_vector& to_sub)
  {
    if(to_sub.size()==1)
      return *this-=to_sub.first();
    else
      {
        if(size()!=to_sub.size())
          {std::cout<<"Error in MC_double_vector::operator-=(MC_double_vector), size are not compatible "<<size()<<";"<<to_sub.size()<<std::endl;exit(-1);}
        
        int N=size();
        for(int k=0;k<N;++k)
          double_list[k]-=to_sub[k];
        return *this;
      }

  }
  MC_double_vector& MC_double_vector::operator*=(const MC_double_vector& to_mult)
  {
    if(to_mult.size()==1)
      return *this*=to_mult.first();
    else
      {
        if(size()!=to_mult.size())
          {std::cout<<"Error in MC_double_vector::operator*=(MC_double_vector), size are not compatible "<<size()<<";"<<to_mult.size()<<std::endl;exit(-1);}
        
        int N=size();
        for(int k=0;k<N;++k)
          double_list[k]*=to_mult[k];
        return *this;
      }
  }
  MC_double_vector& MC_double_vector::operator/=(const MC_double_vector& to_subdiv)
  {
    if(to_subdiv.size()==1)
      return *this/=to_subdiv.first();
    else
      {
        if(size()!=to_subdiv.size())
          {std::cout<<"Error in MC_double_vector::operator/=(MC_double_vector), size are not compatible "<<size()<<";"<<to_subdiv.size()<<std::endl;exit(-1);}
        
        int N=size();
        for(int k=0;k<N;++k)
          {
            double current=to_subdiv[k];
            double epsilon=0.0000001;
            if(fabs(current)<epsilon)
              {std::cout<<"Error in MC_double_vector::operator/=(MC_double_vector), try do divide by u["<<k<<"]=0"<<std::endl;exit(-1);}
            double_list[k]*=current;
          }
        return *this;
      }
  }

  const double& MC_double_vector::first() const
  {
    if(size()>0)
      return double_list[0];
    else
      {
        std::cout<<"Error, call MC_double_vector::first() with size=0"<<std::endl;exit(-1);
        exit(-1);
      }
  }
  double& MC_double_vector::first()
  {
    if(size()>0)
      return double_list[0];
    else
      {
        std::cout<<"Error, call MC_double_vector::first() with size=0"<<std::endl;exit(-1);
        exit(-1);
      }
  }
  const double& MC_double_vector::last() const
  {
    if(size()>0)
      return double_list[size()-1];
    else
      {
        std::cout<<"Error, call MC_double_vector::last() with size=0"<<std::endl;exit(-1);
      }
  }
  double& MC_double_vector::last()
  {
    if(size()>0)
      return double_list[size()-1];
    else
      {
        std::cout<<"Error, call MC_double_vector::last() with size=0"<<std::endl;exit(-1);
      }
  }
  
  std::pair <MC_double_vector,double> MC_double_vector::delete_index(const int& index_to_delete) const
  {
    if(index_to_delete<0 || index_to_delete>=size())
      {std::cout<<"Error in MC_double_vector::delete("<<index_to_delete<<", when size="<<size()<<std::endl;exit(-1);}
    MC_double_vector new_vec;
    int N=size();
    for(int k=0;k<N;++k)
      if(k!=index_to_delete)
        new_vec.add(double_list[k]);
    return std::pair <MC_double_vector,double>(new_vec,double_list[index_to_delete]);
  }
  std::pair <MC_double_vector,MC_double_vector> MC_double_vector::delete_index(const MC_int_vector& index_to_delete) const
  {
    std::set <int> set_index=index_to_delete.to_set();
    MC_double_vector new_vec;
    MC_double_vector deleted_vec;
    int N=size();
    for(int k=0;k<N;k++)
      {
        if(set_index.find(k)!=set_index.end())
          new_vec.add(double_list[k]);
        else
          deleted_vec.add(double_list[k]);
      }
    return std::pair <MC_double_vector,MC_double_vector> (new_vec,deleted_vec);
  }
    
  MC_double_vector MC_double_vector::operator()(const MC_int_vector& index) const
  {
    int N=size();
    int N2=index.size();
    MC_double_vector new_vec=zeros(N2);
    for(int k=0;k<N2;++k)
      {
        if(index[k]<0 || index[k]>N)
          {std::cout<<"Error in MC_double_vector::operator()(MC_double_vector), value index["<<k<<"]="<<index[k]<<" not correct for size="<<N<<std::endl;exit(-1);}
        new_vec[k]=double_list[index[k]];
      }
    return new_vec;
  }
  MC_double_vector MC_double_vector::operator[](const MC_int_vector& index) const
  {return (*this)(index);}


  MC_double_vector& operator>>(std::istream& input,MC_double_vector& vec)
  {
      while(input.good())
      {
          double temp=0;
          input>>temp;
          if(input.good())
              vec.add(temp);
      }
      return vec;
  }


  MC_int_vector operator<(const MC_double_vector& vec,const double& val)
  {
      int N=vec.size();
      MC_int_vector res=MC_int_vector::zeros(N);
      for(int k=0;k<N;k++)
          if(vec[k]<val)
              res[k]=1;
      return res;
  }
  MC_int_vector operator<=(const MC_double_vector& vec,const double& val)
  {
      int N=vec.size();
      MC_int_vector res=MC_int_vector::zeros(N);
      for(int k=0;k<N;k++)
          if(vec[k]<=val)
              res[k]=1;
      return res;
  }
  MC_int_vector operator>(const MC_double_vector& vec,const double& val)
  {
      int N=vec.size();
      MC_int_vector res=MC_int_vector::zeros(N);
      for(int k=0;k<N;k++)
          if(vec[k]>val)
              res[k]=1;
      return res;
  }
  MC_int_vector operator>=(const MC_double_vector& vec,const double& val)
  {
      int N=vec.size();
      MC_int_vector res=MC_int_vector::zeros(N);
      for(int k=0;k<N;k++)
          if(vec[k]>=val)
              res[k]=1;
      return res;
  }
  MC_int_vector operator==(const MC_double_vector& vec,const double& val)
  {
      int N=vec.size();
      MC_int_vector res=MC_int_vector::zeros(N);
      double epsilon=0.0000001;
      for(int k=0;k<N;k++)
          if(fabs(vec[k]-val)<epsilon)
              res[k]=1;
      return res;
  }

  MC_int_vector operator<(const double& val,const MC_double_vector& vec)
  {return vec>val;}
  MC_int_vector operator<=(const double& val,const MC_double_vector& vec)
  {return vec>=val;}
  MC_int_vector operator>(const double& val,const MC_double_vector& vec)
  {return vec<val;}
  MC_int_vector operator>=(const double& val,const MC_double_vector& vec)
  {return vec<=val;}
  MC_int_vector operator==(const double& val,const MC_double_vector& vec)
  {return vec==val;}

  const double* MC_double_vector::pointer() const
  {return (&double_list[0]);}

   MC_double_vector MC_double_vector::linspace(const double& begin,const double& end,const double& iteration)
   {
       if( (begin<end && iteration<0) || (begin>end && iteration>0) || fabs(iteration)<0.0000001)
       {std::cout<<"Error in MC_double_vector::linspace("<<begin<<","<<end<<","<<iteration<<"), loop will never end"<<std::endl;exit(-1);}

       MC_double_vector temp;
       for(double u=begin;u<=(end+iteration/5.0);u+=iteration)
           temp.add(u);
       if(temp.size()>2)
           temp[temp.size()-1]=end; //avoid floating point problem
       return temp;
   }
   MC_double_vector MC_double_vector::sample(const double& begin,const double& end,const int& N_sample)
   {
       if(N_sample<1){std::cout<<"Error in MC_double_vector::sample("<<begin<<","<<end<<","<<N_sample<<"), N_sample must be>=1"<<std::endl;exit(-1);}
       if(N_sample==1)
           return begin;

       double iteration=(end-begin)/static_cast<double>(N_sample-1);
       MC_double_vector res=MC_double_vector::zeros(N_sample-1);
       res.add(begin);
       for(int k=1;k<N_sample;++k)
           res[k]=res[k-1]+iteration;
       res(N_sample-1)=end;
       return res;
   }

   MC_v3d_vector operator*(const MC_double_vector& to_mult,const MC_v3d& vec)
   {
       int N=to_mult.size();
       MC_v3d_vector res(N);
       for(int k=0;k<N;++k)
           res[k]=vec*to_mult[k];
       return res;
   }
   MC_v3d_vector operator*(const MC_v3d& to_mult,const MC_double_vector& vec)
   {return vec*to_mult;}
   MC_v3d_vector operator/(const MC_v3d& vec,const MC_double_vector& to_subdiv)
   {
       double epsilon=0.000001;

       int N=to_subdiv.size();
       MC_v3d_vector res(N);
       for(int k=0;k<N;++k)
       {
           if(fabs(to_subdiv[k])>epsilon)
               res[k]=vec/to_subdiv[k];
           else
           {std::cout<<"Error in MC_v3d_vector::operator/("<<vec<<" MC_double_vector), at position "<<k<<", value MC_double_vector["<<k<<"] is zero"<<std::endl;exit(-1);}
       }
       return res;
   }

    std::pair <double,int> MC_double_vector::min() const
    {

        if(size()<=0)
        {std::cout<<"Error in MC_double_vector::min(), size()="<<size()<<std::endl;exit(-1);}

        double current_min=first();
        int current_index=0;

        int N=size();
        for(int k=0;k<N;++k)
            if(current_min>double_list[k])
            {current_min=double_list[k];current_index=k;}

        return std::pair <double,int> (current_min,current_index);
    }
    std::pair <double,int> MC_double_vector::max() const
    {
        if(size()<=0)
        {std::cout<<"Error in MC_double_vector::max(), size()="<<size()<<std::endl;exit(-1);}

        double current_max=first();
        int current_index=0;

        int N=size();
        for(int k=0;k<N;++k)
            if(current_max<double_list[k])
            {current_max=double_list[k];current_index=k;}

        return std::pair <double,int> (current_max,current_index);
    }

    double MC_double_vector::sum(const MC_double_vector& input)
    {
        double s=0.0;
        int N=input.size();
        for(int k=0;k<N;++k)
            s+=input[k];
        return s;
    }

    std::pair <MC_double_vector,MC_int_vector> MC_double_vector::sort() const
    {
        //prepare the sorting structure
        int N=size();
        std::vector <std::pair<double,int> > vec_to_sort(N);
        for(int k=0;k<N;++k)
            vec_to_sort[k]=std::pair <double,int>(double_list[k],k);

        //sort
        std::sort(vec_to_sort.begin(),
                  vec_to_sort.end(),
                  comparator_less_double_int());

        //redistribute
        std::pair <MC_double_vector,MC_int_vector> res;
        res.first.resize(N);res.second.resize(N);
        for(int k=0;k<N;++k)
        {
            res.first[k]=vec_to_sort[k].first;
            res.second[k]=vec_to_sort[k].second;
        }
        return res;
    }

    MC_double_vector::MC_double_vector(const MC_v3d& vec)
    {
        double_list.push_back(vec[0]);
        double_list.push_back(vec[1]);
        double_list.push_back(vec[2]);
    }
    double* MC_double_vector::pointer_unprotected()
    {return &double_list[0];}

    MC_double_vector MC_double_vector::abs(const MC_double_vector& input)
    {
        int N=input.size();
        MC_double_vector res=MC_double_vector::zeros(N);
        for(int k=0;k<N;++k)
            res[k]=std::abs(input[k]);
        return res;
    }

    MC_double_vector operator*(const double& to_mult,const MC_double_vector& vec)
    {return vec*to_mult;}

//    MC_double_vector MC_double_vector::pointwise_mult(const MC_double_vector& p) const
//    {
//        if(p.size()!=size())
//        {std::cout<<"Error in MC_double_vector::pointwise_mult(), size is not correct:"<<p.size()<<","<<size()<<std::endl;exit(-1);}
//
//        MC_double_vector v=MC_double_vector::zeros(p.size());
//        for(int k=0,N=p.size();k<N;++k)
//            v[k]=(*this)[k]*p[k];
//
//        return v;
//    }

    double MC_double_vector::norm_2() const
    {
        double n=0.0;
        for(int k=0,N=size();k<N;++k)
        {
            double a=double_list[k];
            n += a*a;
        }
        return std::sqrt(n);
    }
    double MC_double_vector::dot(const MC_double_vector& v) const
    {
        int N=size();
        if(N!=v.size())
        {std::cout<<"Error in MC_double_vector::dot(const MC_double_vector& v), size are not consistent: "<<N<<" "<<v.size()<<std::endl;exit(-1);}

        const double* p_v=v.pointer();

        double res=0.0;
        for(int k=0;k<N;++k)
            res += double_list[k]*p_v[k];

        return res;
    }

    MC_double_vector MC_double_vector::inverted() const
    {
        unsigned int N=size();
        MC_double_vector inv=MC_double_vector::zeros(N);
        for(int k=0;k<N;++k)
            inv[k]=double_list[N-k-1];
        return inv;
    }

}