// GenFem - A high-level finite element library
// Copyright (C) 2010-2026 Eric Bechet
//
// See the LICENSE file for license information and contributions.
// Please report all bugs and problems to <bechet@cadxfem.org>.
//
// Initial design: Eric Bechet (rev.734)


#include "genSimpleFunctionExtended.h"
#include "GModel.h"
#include "genDataIO.h"

#include <iostream>
#include <iomanip>
#include <fstream>
#include <cassert>
#include <boost/algorithm/string.hpp>

//const char genBoundaryCondition::physNames[][256]={"Displacement","Force","Potential","Charge","Temperature","ThermalFlux"};

genBoundaryCondition::~genBoundaryCondition()
{
  if (own)
  {
    if (fscalar) delete fscalar;
    if (fvector) delete fvector;
    if (ftensor) delete ftensor;
    if (fcscalar) delete fcscalar;
    if (fcvector) delete fcvector;
    if (fctensor) delete fctensor;
    if (periodic) delete periodic;
  }
}


genDataIO::~genDataIO()
{
  for (int i=0;i<Domains.size();++i) delete Domains[i];
  for (int i=0;i<BCs.size();++i) delete BCs[i];
  if (Model) delete Model;
}

void genDataIO::readInputFile(const std::string &fileName)
{
  std::ifstream ifs(fileName.c_str());
  if (!ifs)
  {
    printf("Unable to open file '%s'\n", fileName.c_str());
    exit(0);
  }
  while ( !ifs.eof() )
  {
    std::string token;
    NextTokenWithoutComments(ifs, token);
    readToken(ifs, token);
  }
  ifs.close();
}

void genDataIO::readToken(std::istream &is, std::string &token)
{
  if ( token=="Dim")
  {
    readScalar(is,dim);
  }
  else if (token=="MeshFile")
  {
    readString(is,GModelBaseName);
    std::string::size_type idx = GModelBaseName.find(".msh");
    GModelBaseName = GModelBaseName.substr(0, idx).c_str();

    std::string GModelFileName;
    GModelFileName = GModelBaseName + ".msh";

    printf("mesh name : %s\n", GModelFileName.c_str());
    Model = new GModel();
    int readingDone;
    readingDone = Model->readMSH(GModelFileName.c_str());
    if (!readingDone) exit(0);

    if (dim!=0) dim = dim;
    else if (Model->getNumRegions()) dim = 3;
    else if (Model->getNumFaces())   dim = 2;
    else if (Model->getNumEdges())   dim = 1;
  }
  else if (token=="BoundaryCondition")
  {
    genBoundaryCondition* BC = new genBoundaryCondition;
    is >> *BC;
    BCs.push_back ( BC );
  }
  else if (token=="Domain" )
  {
    genDomain* Sup=new genDomain;
    is >> *Sup;
    Domains.push_back(Sup);
  }
  else if ( token=="Solver")
  {
    readScalar(is,solvertype);
  }
  else if (token.find_first_of("/#.;,\\+-*[]")==std::string::npos)
  {
    std::string tenstype;
    NextTokenWithoutComments(is,tenstype);
    if (tenstype=="Integer")
    {
      readScalar(is,User.Integers[token]);
      User.AllData[token]=genDomain::Integer;
    }
    else if (tenstype=="Scalar")
    {
      readScalar(is,User.Scalars[token]);
      User.AllData[token]=genDomain::Scalar;
    }
    else if (tenstype=="CScalar")
    {
      readScalar(is,User.CScalars[token]);
      User.AllData[token]=genDomain::CScalar;
    }
    else if (tenstype=="Vector")
    {
      readVector(is,User.Vectors[token]);
      User.AllData[token]=genDomain::Vector;
    }
    else if (tenstype=="CVector")
    {
      readVector(is,User.CVectors[token]);
      User.AllData[token]=genDomain::CVector;
    }
    else if (tenstype=="Matrix")
    {
      readMatrix(is,User.Tensors[token]);
      User.AllData[token]=genDomain::Tensor;
    }
    else if (tenstype=="CMatrix")
    {
      readMatrix(is,User.CTensors[token]);
      User.AllData[token]=genDomain::CTensor;
    }

  }
}



void genDataIO::dumpContent(std::ostream &os)
{
  if (Model)
  {
    os << "MeshFile " << GModelBaseName << std::endl;
  }
  for (int i=0;i<Domains.size();++i) os << *Domains[i] << std::endl;
  for (int i=0;i<BCs.size();++i) os << *BCs[i] << std::endl;
  for (std::map<std::string,genDomain::tenstype>::const_iterator it=User.AllData.begin();it!=User.AllData.end();++it)
  {
    if (it->second==genDomain::Integer)
      os << it->first << " Integer " << User.Integers.find(it->first)->second << std::endl;
    else if (it->second==genDomain::Scalar)
      os << it->first << " Scalar " << User.Scalars.find(it->first)->second << std::endl;
    else if (it->second==genDomain::CScalar)
      os << it->first << " CScalar " << User.CScalars.find(it->first)->second << std::endl;
    else if (it->second==genDomain::Vector)
      os << it->first << " Vector " << "XXXXXXXXX"/*User.Vectors.find(it->first)->second */<< std::endl;
    else if (it->second==genDomain::CVector)
      os << it->first << " CVector " << "XXXXXXXXX"/*User.CVectors.find(it->first)->second */<< std::endl;
    else if (it->second==genDomain::Tensor)
      os << it->first << " Matrix " << "XXXXXXXXX"/*User.Tensors.find(it->first)->second */<< std::endl;
    else if (it->second==genDomain::CTensor)
      os << it->first << " CMatrix " << "XXXXXXXXX"/*User.CTensors.find(it->first)->second */<< std::endl;
  }
  os << "SolverType " << solvertype << std::endl;
}


std::istream & operator>>(std::istream &is, genDomain &obj)
{
  int physical;
  int dim;
  readScalar(is,physical);
  readScalar(is,dim);
  obj.tag = physical;
  obj.dim = dim;
  NextTokenWithoutComments(is,obj.type);
  std::string token;
  NextTokenWithoutComments(is,token);
  if (token=="Data")
  {
    std::string name;
    NextTokenWithoutComments(is,name);
    while (name!="EndData")
    {
      std::string tenstype;
      NextTokenWithoutComments(is,tenstype);
      if (tenstype=="Integer")
      {
        readScalar(is,obj.Integers[name]);
        obj.AllData[name]=genDomain::Integer;
      } 
      else if (tenstype=="Scalar")
      {
        readScalar(is,obj.Scalars[name]);
        obj.AllData[name]=genDomain::Scalar;
      }
      else if (tenstype=="CScalar")
      {
        readScalar(is,obj.CScalars[name]);
        obj.AllData[name]=genDomain::CScalar;
      }
      else if (tenstype=="Vector")
      {
        readVector(is,obj.Vectors[name]);
        obj.AllData[name]=genDomain::Vector;
      }
      else if (tenstype=="CVector")
      {
        readVector(is,obj.CVectors[name]);
        obj.AllData[name]=genDomain::CVector;
      }
      else if (tenstype=="Matrix")
      {
        readMatrix(is,obj.Tensors[name]);
        obj.AllData[name]=genDomain::Tensor;
      }
      else if (tenstype=="CMatrix")
      {
        readMatrix(is,obj.CTensors[name]);
        obj.AllData[name]=genDomain::CTensor;
      }
      NextTokenWithoutComments(is,name);
    }
  }
  return is;
}

std::ostream & operator<<(std::ostream &os, const genDomain &obj)
{
  os << "Domain " << obj.tag << " " << obj.dim << " " << obj.type << std::endl;
  os << "Data" << std::endl;
  for (std::map<std::string,genDomain::tenstype>::const_iterator it=obj.AllData.begin();it!=obj.AllData.end();++it)
  {
    if (it->second==genDomain::Integer)
      os << it->first << " Integer " << obj.Integers.find(it->first)->second << std::endl;
    else if (it->second==genDomain::Scalar)
      os << it->first << " Scalar " << obj.Scalars.find(it->first)->second << std::endl;
    else if (it->second==genDomain::CScalar)
      os << it->first << " CScalar " << obj.CScalars.find(it->first)->second << std::endl;
    else if (it->second==genDomain::Vector)
      os << it->first << " Vector " << "XXXXXXXXX"/*obj.Vectors.find(it->first)->second */<< std::endl;
    else if (it->second==genDomain::CVector)
      os << it->first << " CVector " << "XXXXXXXXX"/*obj.Vectors.find(it->first)->second */<< std::endl;
    else if (it->second==genDomain::Tensor)
      os << it->first << " Matrix " << "XXXXXXXXX"/*obj.Tensors.find(it->first)->second */<< std::endl;
    else if (it->second==genDomain::CTensor)
      os << it->first << " CMatrix " << "XXXXXXXXX"/*obj.Tensors.find(it->first)->second */<< std::endl;
  }
  os << "EndData";
  return os;
}


std::istream & operator>>(std::istream &is, genBoundaryCondition &obj)
{
  std::string topo;
  obj.dim=genBoundaryCondition::UNDEF;
  obj.kind="";
  do
  {
    NextTokenWithoutComments(is,topo);
    if (topo=="Node") obj.dim=genBoundaryCondition::ON_VERTEX;
    else if (topo=="Edge") obj.dim=genBoundaryCondition::ON_EDGE;
    else if (topo=="Face") obj.dim=genBoundaryCondition::ON_FACE;
    else if (topo=="Volume") obj.dim=genBoundaryCondition::ON_VOLUME;
    else obj.kind=topo;
  } while ((obj.dim==genBoundaryCondition::UNDEF)&&!is.eof());
  readScalar(is,obj.tag);
  std::string field;
//  obj.What =genBoundaryCondition::UserDefined;
  NextTokenWithoutComments(is,obj.What); // DISPLACEMENT , POTENTIAL, etc...
//  for (int i=0;i<(int)genBoundaryCondition::Periodic;++i)
//    if (obj.user == genBoundaryCondition::physNames[i] ) {obj.What=(genBoundaryCondition::physical)i;break;}
  std::string nnn;
  std::string iswhat;
  bool ok=false;
  do
  {
    NextTokenWithoutComments(is,nnn);
    if (nnn=="Constant")
    {
      NextTokenWithoutComments(is,iswhat); // Scalar,Vector
      if (iswhat=="Scalar")
      {
        obj.isWhat=genBoundaryCondition::Scalar;
        double val;
        readScalar(is,val);
        obj.fscalar = new genSimpleFunctionConstant<genTensor0<double> > ( val );
      }
      if (iswhat=="CScalar")
      {
        obj.isWhat=genBoundaryCondition::CScalar;
        std::complex<double> val;
        readScalar(is,val);
        obj.fcscalar = new genSimpleFunctionConstant<genTensor0<std::complex<double> > > ( val );
      }
      else if (iswhat=="Vector")
      {
        obj.isWhat=genBoundaryCondition::Vector;
        std::vector<double> val;
        readVector(is,val);
        obj.fvector = new genSimpleFunctionConstant<genTensor1<double> > ( genTensor1<double> ( val[0], val[1], val[2] ) );
      }
      else if (iswhat=="CVector")
      {
        obj.isWhat=genBoundaryCondition::CVector;
        std::vector<std::complex<double> > val;
        readVector(is,val);
        obj.fcvector = new genSimpleFunctionConstant<genTensor1<std::complex<double> > > ( genTensor1<std::complex<double> > ( val[0], val[1], val[2] ) );
      }
      else if (iswhat=="Matrix")
      {
        obj.isWhat=genBoundaryCondition::Tensor;
        genMatrix<double> val;
        readMatrix(is,val);
        genTensor2<double> M; M.setMat(val);
        obj.ftensor = new genSimpleFunctionConstant<genTensor2<double> > ( M );
      }
      else if (iswhat=="CMatrix")
      {
        obj.isWhat=genBoundaryCondition::CTensor;
        genMatrix<std::complex<double> > val;
        readMatrix(is,val);
        genTensor2<std::complex<double> > M; M.setMat(val);
        obj.fctensor = new genSimpleFunctionConstant<genTensor2<std::complex<double> > > ( M );
      }
      ok=false;
    }
    else
    if (nnn=="Function")
    {
      NextTokenWithoutComments(is,iswhat); // Scalar,Vector
      if (iswhat=="Scalar")
      {
          obj.isWhat=genBoundaryCondition::Scalar;
          std::vector<std::string> expr(1),var;
          readString(is,expr[0]);
          readVector(is,var);
          obj.fscalar = new genSimpleFunctionAnalytical<genTensor0<double> >(expr,var);
      }
      if (iswhat=="CScalar")
      {
          obj.isWhat=genBoundaryCondition::CScalar;
          std::vector<std::string> expr(1),var;
          readString(is,expr[0]);
          readVector(is,var);
//          obj.fcscalar = new genSimpleFunctionAnalytical<genTensor0<std::complex<double> > >(expr,var);
      }
      else if (iswhat=="Vector")
      {
        obj.isWhat=genBoundaryCondition::Vector;
        std::vector<std::string> expr,var;
        readVector(is,expr);
        readVector(is,var);
        obj.fvector = new genSimpleFunctionAnalytical<genTensor1<double> >(expr,var);
      }
      else if (iswhat=="CVector")
      {
        obj.isWhat=genBoundaryCondition::CVector;
        std::vector<std::string> expr,var;
        readVector(is,expr);
        readVector(is,var);
//        obj.fcvector = new genSimpleFunctionAnalytical<genTensor1<std::complex<double> > >(expr,var);
      }
      ok=false;
    }
    else
    if (nnn=="Periodic")
    {
      obj.kind="Periodic";
      obj.periodic=new genBoundaryCondition;
      is>>*(obj.periodic);
      std::vector<std::string> expr,var;
      readVector(is,expr);
      readVector(is,var);
      obj.fvector = new genSimpleFunctionAnalytical<genTensor1<double> >(expr,var);
      ok=false;
    }
    else
    {
      obj.comp1=atoi(nnn.c_str());
      ok=true;
    }
  } while (ok);
  return is;
}

std::ostream & operator<<(std::ostream &os, const genBoundaryCondition &obj)
{
  os << "BoundaryCondition ";
  if (obj.dim==genBoundaryCondition::ON_VERTEX) os << "Vertex ";
  else if (obj.dim==genBoundaryCondition::ON_EDGE) os << "Edge ";
  else if (obj.dim==genBoundaryCondition::ON_FACE) os << "Face ";
  else if (obj.dim==genBoundaryCondition::ON_VOLUME) os << "Volume ";
  os << obj.tag << " " << obj.isWhat ;
  if (obj.periodic!=NULL) os << " Periodic ";

  return os;
}



void NextToken(std::istream &is, std::string &token, const std::string &next)
{
  if (next=="")
    is >> token;
  else
    token=next;
}
void RemoveComments(std::istream &is, std::string &token)
{
  size_t found;
  found = token.find("//");
  if (found != std::string::npos)
  {
    std::string comment;
    std::getline(is, comment);
//     comment = token.substr(found) + comment;
    token = token.substr(0,found);
//     printf("  comment=\"%s\"\n", comment.c_str());
  }
}
void NextTokenWithoutComments(std::istream &is, std::string &token, const std::string &next)
{
  NextToken(is, token, next);
  RemoveComments(is, token);
  if ( (token=="") && (!is.eof()) )
    NextTokenWithoutComments(is, token);
//   else if (!is.eof())
//     printf("token=\"%s\"\n", token.c_str());
}


void readScalar(std::istream &is, int &s)
{
  std::string str;
  NextTokenWithoutComments(is, str);
  s = atoi(str.c_str());
}

void readScalar(std::istream &is, double &s)
{
  std::string str;
  NextTokenWithoutComments(is, str);
  s = atof(str.c_str());
}


void readString(std::istream &is, std::string &s)
{
  NextTokenWithoutComments(is, s);
}

void readArrayString(std::istream &is, std::string &s)
{
  s.clear();

  std::string token;
  NextTokenWithoutComments(is, token);

  size_t foundBegin = token.find("[");
  size_t foundEnd = token.find("]");

  if (foundBegin == 0)
  {
    while ( (foundEnd == std::string::npos) && (!is.eof()) )
    {
      s += token; 
      NextTokenWithoutComments(is, token);
      foundEnd = token.find("]");
    }
    s += token;
  }

  if ( (foundBegin != 0) ||  (foundEnd != token.size()-1) || (is.eof()) )
    printf("Error in readArrayString : \"%s\"\n", s.c_str());
//   else
//     printf("Result of readArrayString : \"%s\"\n", s.c_str());
}


void readParArrayString(std::istream &is, std::string &s)
{
  s.clear();

  std::string token;
  NextTokenWithoutComments(is, token);

  size_t foundBegin = token.find("(");
  size_t foundEnd = token.find(")");

  if (foundBegin == 0)
  {
    while ( (foundEnd == std::string::npos) && (!is.eof()) )
    {
      s += token; 
      NextTokenWithoutComments(is, token);
      foundEnd = token.find(")");
    }
    s += token;
  }

  if ( (foundBegin != 0) ||  (foundEnd != token.size()-1) )
    s+="("+token+", 0.0)"; // single real value
    
    if (is.eof())
    printf("Error in readParArrayString : \"%s\"\n", s.c_str());
//   else
//     printf("Result of readParArrayString : \"%s\"\n", s.c_str());
}

void readParVector(const std::string &string_v, std::vector<double> &v)
{
  assert(string_v.find('(')==0);
  assert(string_v.find(')')==string_v.size()-1);

  std::string s = string_v.substr(1,string_v.size()-2);
  std::vector<std::string> comp;

  boost::split( comp, s, boost::is_any_of(",") );

  assert(comp.size()!=0);

  if (v.size() != comp.size() )
    v.resize(comp.size());
  for (int i=0; i<comp.size(); ++i)
    v[i] = atof(comp[i].c_str());
}

void readScalar(std::istream &is, std::complex<double> &s)
{
  std::string str;
  std::vector<double> v;
  readParArrayString(is,str);
  readParVector(str,v);
  if (v.size()) s.real(v[0]); else s.real(0.);
  if (v.size()>1) s.imag(v[1]); else s.imag(0.);
}

void readVector(std::istream &is, std::vector<int> &v)
{
  std::string s;
  readArrayString(is,s);
  readVector(s,v);
}

void readVector(std::istream &is, std::vector<double> &v)
{
  std::string s;
  readArrayString(is,s);
  readVector(s,v);
}

void readVector(std::istream &is, std::vector<std::complex<double> > &v)
{
  std::string s;
  readArrayString(is,s);
  readVector(s,v);
}

void readVector(std::istream &is, std::vector<std::string> &v)
{
  std::string s;
  readArrayString(is,s);
  readVector(s,v);
}


void readMatrix(std::istream &is, genMatrix<int> &m)
{
  std::string s;
  readArrayString(is,s);
  readMatrix(s,m);
}

void readMatrix(std::istream &is, genMatrix<double> &m)
{
  std::string s;
  readArrayString(is,s);
  readMatrix(s,m);
}

void readMatrix(std::istream &is, genMatrix<std::complex<double> > &m)
{
  std::string s;
  readArrayString(is,s);
  readMatrix(s,m);
}

void readMatrix(std::istream &is, genMatrix<std::string> &m)
{
  std::string s;
  readArrayString(is,s);
  readMatrix(s,m);
}

void readVector(const std::string &string_v, std::vector<int> &v)
{
  assert(string_v.find('[')==0);
  assert(string_v.find(']')==string_v.size()-1);

  std::string s = string_v.substr(1,string_v.size()-2);
  std::vector<std::string> comp;

  boost::split( comp, s, boost::is_any_of(",") );

  assert(comp.size()!=0);

  if (v.size() != comp.size() )
    v.resize(comp.size());
  for (int i=0; i<comp.size(); ++i)
    v[i] = atof(comp[i].c_str());
}

void readVector(const std::string &string_v, std::vector<double> &v)
{
  assert(string_v.find('[')==0);
  assert(string_v.find(']')==string_v.size()-1);

  std::string s = string_v.substr(1,string_v.size()-2);
  std::vector<std::string> comp;

  boost::split( comp, s, boost::is_any_of(",") );

  assert(comp.size()!=0);

  if (v.size() != comp.size() )
    v.resize(comp.size());
  for (int i=0; i<comp.size(); ++i)
    v[i] = atof(comp[i].c_str());
}

void readComplex(const std::string &comp,std::complex<double> &val)
{
  if (comp.find('(')!=std::string::npos) //complex
  {
    std::vector<std::string> s;
    boost::split(s, comp, boost::is_any_of("(,)") );
//     std::cout << comp << "!" << std::endl;
//     std::cout << s.size()<< std::endl;
    bool first=true;
    val.real(0.);
    val.imag(0.);
    for (int i=0;i<s.size();++i)
    {
//       std::cout << s[i] << std::endl;

      if (s[i].size()!=0)
      {
        if (first)
        {
          val.real(atof(s[i].c_str()));
          first=false;
        }
        else
        {
          val.imag(atof(s[i].c_str()));
          break;
        }
      }
    }
  }
  else // real
  {
    val.real(atof(comp.c_str())) ;
    val.imag(0.0);
  }
}

void readVector(const std::string &string_v, std::vector<std::complex<double> > &v)
{
  assert(string_v.find('[')==0);
  assert(string_v.find(']')==string_v.size()-1);

  std::string s = string_v.substr(1,string_v.size()-2);
  std::string comp;
  v.clear();
  std::complex<double> val;
  int i=0;
  int ii=0;
  bool c;
  while (i<s.size())
  {
    ii=i;
    while (s[i]!=',')
    {
      if (i==s.size()) break;
      if (s[i]=='(' )
      {
        while (s[i]!=')') {i++;if (i==s.size()) break;}
      }
      i++;
    }
    comp=s.substr(ii,i-ii);
//    std::cout << comp << std::endl;
    readComplex(comp,val);
    v.push_back(val);
    comp.clear();
    i++;
  }
//  std::cout << std::endl;
//  for (int i=0;i<v.size();++i) std::cout << v[i] << " " ;
//  std::cout << std::endl;
}

void readVector(const std::string &string_v, std::vector<std::string> &v)
{
  assert(string_v.find('[')==0);
  assert(string_v.find(']')==string_v.size()-1);

  std::string s = string_v.substr(1,string_v.size()-2);
  std::vector<std::string> comp;

  boost::split( comp, s, boost::is_any_of(",") );

  assert(comp.size()!=0);

  v = comp;
}


void readArray(const std::string &string_v, double t[], int size)
{
  std::vector<double> v;
  readVector(string_v, v);

  assert(v.size()==size);
  std::copy( v.begin(), v.end(), t);
}

void readArray(const std::string &string_v, std::complex<double> t[], int size)
{
  std::vector<std::complex<double> > v;
  readVector(string_v, v);

  assert(v.size()==size);
  std::copy( v.begin(), v.end(), t);
}

void readMatrix(const std::string &string_m, genMatrix<int> &m)
{
  assert(string_m.find('[')==0);
  assert(string_m.find(']')==string_m.size()-1);

  std::string s =string_m.substr(1,string_m.size()-2);
  std::vector<std::string> string_line;
  std::vector< std::vector<std::string> > comp;
  int size_i, size_j;

  boost::split( string_line, s, boost::is_any_of(";") );
  comp.resize(string_line.size());
  size_i=string_line.size();

  for (int i=0; i<size_i; ++i)
  {
    std::vector<std::string> temp;
    boost::split( temp, string_line[i], boost::is_any_of(",") );
    comp[i].resize(temp.size());
    for (int j=0; j<temp.size(); ++j)
      comp[i][j]=temp[j];

    if (i == 0) size_j=comp[i].size();
    else assert(size_j == comp[i].size());
  }
  m.resize(size_i, size_j);

  for (int i=0; i<size_i; ++i)
  {
    for (int j=0; j<size_j; ++j)
    {
      m(i,j) = atoi(comp[i][j].c_str());
    }
  }
}

void readMatrix(const std::string &string_m, genMatrix<double> &m)
{
  assert(string_m.find('[')==0);
  assert(string_m.find(']')==string_m.size()-1);
  std::string s =string_m.substr(1,string_m.size()-2);
  std::vector<std::string> string_line;
  std::vector< std::vector<std::string> > comp;
  int size_i, size_j;

  boost::split( string_line, s, boost::is_any_of(";") );
  comp.resize(string_line.size());
  size_i=string_line.size();

  for (int i=0; i<size_i; ++i)
  {
    std::vector<std::string> temp;
    boost::split( temp, string_line[i], boost::is_any_of(",") );
    comp[i].resize(temp.size());
    for (int j=0; j<temp.size(); ++j)
      comp[i][j]=temp[j];

    if (i == 0) size_j=comp[i].size();
    else assert(size_j == comp[i].size());
  }
  m.resize(size_i, size_j);

  for (int i=0; i<size_i; ++i)
  {
    for (int j=0; j<size_j; ++j)
    {
      m(i,j) = atof(comp[i][j].c_str());
    }
  }
}

void readMatrix(const std::string &string_m, genMatrix<std::complex<double > > &m) // works only for reals as of now.
{
  assert(string_m.find('[')==0);
  assert(string_m.find(']')==string_m.size()-1);
  std::string s =string_m.substr(1,string_m.size()-2);
  std::vector<std::string> string_line;
  std::vector< std::vector<std::string> > comp;
  int size_i, size_j;

  boost::split( string_line, s, boost::is_any_of(";") );
  comp.resize(string_line.size());
  size_i=string_line.size();

  for (int i=0; i<size_i; ++i)
  {
    std::vector<std::string> temp;
    boost::split( temp, string_line[i], boost::is_any_of(",") );
    comp[i].resize(temp.size());
    for (int j=0; j<temp.size(); ++j)
      comp[i][j]=temp[j];

    if (i == 0) size_j=comp[i].size();
    else assert(size_j == comp[i].size());
  }
  m.resize(size_i, size_j);

  for (int i=0; i<size_i; ++i)
  {
    for (int j=0; j<size_j; ++j)
    {
      m(i,j) = atof(comp[i][j].c_str());
    }
  }
}

void readMatrix(const std::string &string_m, genMatrix<std::string> &m)
{
  assert(string_m.find('[')==0);
  assert(string_m.find(']')==string_m.size()-1);

  std::string s =string_m.substr(1,string_m.size()-2);
  std::vector<std::string> string_line;
  std::vector< std::vector<std::string> > comp;
  int size_i, size_j;

  boost::split( string_line, s, boost::is_any_of(";") );
  comp.resize(string_line.size());
  size_i=string_line.size();

  for (int i=0; i<size_i; ++i)
  {
    std::vector<std::string> temp;
    boost::split( temp, string_line[i], boost::is_any_of(",") );
    comp[i].resize(temp.size());
    for (int j=0; j<temp.size(); ++j)
      comp[i][j]=temp[j];

    if (i == 0) size_j=comp[i].size();
    else assert(size_j == comp[i].size());
  }

//  std::string* empty = new std::string[size_i*size_j];
//  for (int i=0; i<size_i*size_j; ++i) empty[i]=" ";
  genMatrix<std::string> temp(size_i, size_j);

  for (int i=0; i<size_i; ++i)
  {
    for (int j=0; j<size_j; ++j)
    {
      temp(i,j) = comp[i][j];
    }
  }
  
  m = temp;
}