class NOperators;

class ClusterData{
  function2D<int> F_i;
  function2D<function2D<double> > F_M;
public:
  function2D<double> Ene;
private:
  function2D<double> Ene0;
  function2D<double> Spin;
  function1D<int> Ks;
  function1D<int> msize_;
  function2D<double> Patom_;
  function1D<double> Patom;
  vector<function2D<double> > Ma;
  vector<int> gflip_index;
public:
  int Nm, nsize, N_ifl, N_flavors;
  int max_size, Osize;
  
  function1D<double> Sz;
  function1D<int> ifl_dim;
  function1D<int> ifl_from_fl;
  function1D<int> bfl_from_fl;
  function2D<int> fl_from_ifl;
  vector<function2D<int> > tfl_index;
  vector<deque<int> > sign;
  vector<deque<int> > conjg;
  vector<deque<int> > bfl_index;
  Number Zatom;
  function1D<function2D<function2D<double> > > HF_M;
  function1D<double> natom;
  int N_unique_fl;
  function1D<int> fl_deg;
  function1D<int> Ns;
  function1D<double> epsk;
  string RealSigma, ImagSigma;
  map<int,double> Id;
  function2D<int> gflip_fl;
  function1D<pair<int,int> > gflip_ifl;
  function2D<bool> ifl_equivalent;
  function1D<int> gs_ind;
  bool QHB1;
  int cnsize;
  function1D<int> cNs;
  function1D<int> cmsize_;
  function2D<int> cF_i;
  function1D<int> cindx;
  vector<vector<double> > cEne;
  function2D<function2D<double> > cF_M;
public:
  void ParsInput(ifstream& gout);
  void EvaluateMa(double beta, double mu, double U);
  int size() const {return nsize;}
  int msize(int i) const {return msize_[i];}
  int praState(int i) const {return i+1;}
  int Fi(int iop, int ist) const {return F_i(iop,ist);}
  const funProxy<int>& Fi(int iop) const {return F_i[iop];}
  const function2D<double>& FM(int op, int st) const {return F_M(op,st);}
  const funProxy<function2D<double> >& FM(int op) const {return F_M[op];}
  double patom(int i) const {return Patom[i];}
  int N_baths(int ifl) const {return bfl_index[ifl].size();}
  void RecomputeCix(function2D<double>& AProb, int Naver, double treshold, function1D<NState>& praStates);
  friend class NOperators;
  void Read_for_HB1(ifstream& gout, double mu, double U);
  void HB1(double beta, double mu, const mesh1D& iom, const function2D<double>& AProb, int nom, int aom, const function2D<dcomplex>& Delta, function2D<dcomplex>& Sigma, int nom_small, double sderiv);
  double P_atom(int i, int m) const {return Patom_(i+1,m);}
};

void ClusterData::ParsInput(ifstream& gout)
{
  gout.ignore(1000,'\n');
  gout.ignore(1000,'\n');
  gout>>Nm>>nsize>>N_ifl>>max_size;

  if (!gout || Nm<0 || nsize<0 || N_ifl<0 || max_size<0){cerr<<"Something wrong reading cix file!"<<endl; exit(1);}
  gout.ignore(1000,'\n');
  gout.ignore(1000,'\n');

  ifl_dim.resize(N_ifl);
  sign.resize(N_ifl);
  conjg.resize(N_ifl);
  tfl_index.resize(N_ifl);
  bfl_index.resize(N_ifl);
  gflip_index.resize(N_ifl);
  
  N_flavors=0;
  int max_dim=0;
  for (int i=0; i<N_ifl; i++){
    int it;
    gout>>it>>ifl_dim[i];
    if (it!=i) {cerr<<"Something wrong in reading symmetries of baths in cix file!"<<endl; exit(1);}
    tfl_index[i].resize(ifl_dim[i],ifl_dim[i]);
    int ii=0;
    string str;
    for (int i1=0; i1<ifl_dim[i]; i1++){
      for (int i2=0; i2<ifl_dim[i]; i2++){
	int sgn=1; int cgn=0;
	gout>>str;
	size_t im = str.find("-");
	if (im != string::npos){
	  sgn=-1;
	  str.replace(im,1,"");
	}
	size_t ic = str.find("*");
	if (ic != string::npos){
	  cgn=1;
	  str.replace(ic,1,"");
	}
	int ind = atoi(str.c_str());
	bfl_index[i].push_back(ind);
	sign[i].push_back(sgn);
	conjg[i].push_back(cgn);
	Id[ind] = (i1==i2) ? 1 : 0;
	tfl_index[i][i1][i2] = ii++;
      }
    }
    if (ifl_dim[i]>max_dim) max_dim = ifl_dim[i];
    N_flavors += ifl_dim[i];
    {
      int is;
      gout>>is;
      gflip_index[i]=is;
    }
    gout.ignore(1000,'\n');
  }
  gout.ignore(1000,'\n');

  N_unique_fl = 0;
  for (int ifl=0; ifl<N_ifl; ifl++)
    for (size_t b=0; b<bfl_index[ifl].size(); b++)
      if (bfl_index[ifl][b]>N_unique_fl) N_unique_fl = bfl_index[ifl][b];
  N_unique_fl++;
  
  fl_deg.resize(N_unique_fl);
  fl_deg=0;
  for (int ifl=0; ifl<N_ifl; ifl++)
    for (size_t b=0; b<bfl_index[ifl].size(); b++) fl_deg[bfl_index[ifl][b]]++;

  ifl_from_fl.resize(N_flavors);
  bfl_from_fl.resize(N_flavors);
  fl_from_ifl.resize(N_ifl,max_dim);
  
  int n_flavors=0;
  for (int ifl=0; ifl<N_ifl; ifl++){
    for (int jb=0; jb<ifl_dim[ifl]; jb++){
      ifl_from_fl[n_flavors] = ifl;
      bfl_from_fl[n_flavors] = jb;
      fl_from_ifl(ifl,jb) = n_flavors;
      n_flavors++;
    }
  }

  map<int,deque<int> > gfl_tmp;
  for (int ifl=0; ifl<N_ifl; ifl++){
    int ii = gflip_index[ifl];
    gfl_tmp[ii].push_back(ifl);
  }
  int sz=0;
  for (map<int,deque<int> >::iterator ia=gfl_tmp.begin(); ia!=gfl_tmp.end(); ia++){
    deque<int>& d = ia->second;
    sz += d.size()*(d.size()-1)/2;
  }
  
  gflip_fl.resize(sz,N_flavors);
  gflip_ifl.resize(sz);
  for (int iu=0; iu<sz; iu++){
    for (int j=0; j<gflip_fl.size_Nd(); j++) gflip_fl[iu][j] = j;
  }

  int iu=0;
  for (map<int,deque<int> >::iterator ia=gfl_tmp.begin(); ia!=gfl_tmp.end(); ia++){
    deque<int>& d = ia->second;
    for (size_t j1=0; j1<d.size(); j1++){
      for (size_t j2=j1+1; j2<d.size(); j2++){
	int ifa = d[j1];
	int ifb = d[j2];
	if (ifl_dim[ifa]!=ifl_dim[ifb]){cerr<<"Dimensions of similar baths have to be the same!"<<endl; exit(1);}
	gflip_ifl[iu] = make_pair(ifa,ifb);
	for (int ib=0; ib<ifl_dim[ifa]; ib++){
	  int fa = fl_from_ifl[ifa][ib];
	  int fb = fl_from_ifl[ifb][ib];
	  gflip_fl[iu][fa] = fb;
	  gflip_fl[iu][fb] = fa;
	}
	iu++;
      }
    }
  }
  ifl_equivalent.resize(N_ifl,N_ifl);
  for (int ifl=0; ifl<N_ifl; ifl++){
    for (int jfl=0; jfl<N_ifl; jfl++){
      ifl_equivalent[ifl][jfl] = bfl_index[ifl]==bfl_index[jfl];
    }
  }
  
  if (!gout){cerr<<"Something wrong in reading cix file! Exiting! "<<endl; exit(1);}
  
  //  gout.ignore(1000,'\n');
  epsk.resize(N_unique_fl);
  epsk=0;
  for (int ik=0; ik<N_unique_fl; ik++) gout>>epsk[ik];
  gout.ignore(1000,'\n');
  gout.ignore(1000,'\n');
  
  F_i.resize(2*N_flavors,nsize+1);
  F_M.resize(2*N_flavors,nsize+1);
  Ene.resize(nsize+1,max_size);
  Ene0.resize(nsize+1,max_size);
  Spin.resize(nsize+1,max_size);
  Ns.resize(nsize+1);
  Ks.resize(nsize+1);
  Sz.resize(nsize+1);
  msize_.resize(nsize+1);
  
  msize_=0;
  F_i=0;
  
  int tsize=0;
  for (int i=1; i<=nsize; i++){
    int it1, in, ik, isize;   double dsz;
    gout>>it1>>in>>ik>>dsz>>isize;
    if (it1!=i){cerr<<"Something wrong in parsing cix file!"<<endl;}
    Ns[i] = in;
    Ks[i] = ik;
    Sz[i] = dsz;
    msize_[i] = isize;
    for (int ib=0; ib<N_flavors; ib++){
      int ist;
      gout>>ist;
      if (ist!=0){
	F_i(2*ib,i)=ist;
	F_i(2*ib+1,ist)=i;
      }
    }
    for (int is=0; is<isize; is++) gout>>Ene(i,is);
    Ene[i].resize(isize);
    for (int is=0; is<isize; is++) gout>>Spin(i,is);
    Spin[i].resize(isize);
    gout.ignore(1000,'\n');
    
    tsize += isize;
  }

  gout.ignore(1000,'\n');
  Ene0 = Ene;
  
  for (int i=1; i<=nsize; i++){
    for (int ib=0; ib<N_flavors; ib++){
      int it, jt, size1, size2;
      gout>>it>>jt>>size1>>size2;
      if (it!=i || jt!=F_i(2*ib,i)) cerr<<"Something wrong reading cix file"<<endl;
      F_M(2*ib,i).resize(size2,size1); // constructor changes i->ii
      int ii = F_i(2*ib,i);
      F_M(2*ib+1,ii).resize(size1,size2);// destructor changes ii->i
      for (int i1=0; i1<size1; i1++){
	for (int i2=0; i2<size2; i2++){
	  double m;
	  gout>>m;
	  F_M(2*ib,i)(i2,i1)=m;    // constructor
	  F_M(2*ib+1,ii)(i1,i2)=m; // destructor
	}
      }
      gout.ignore(1000,'\n');
    }
  }
  string str;
  QHB1=false;
  getline(gout,str);
  if (str=="HB1"){
    clog<<"Have HB1"<<endl;
    QHB1=true;
  }else{
    // high frequency expansion, the first few moments of self-energy
    if (gout){
      getline(gout,str);
      RealSigma = str;
      getline(gout,str);
      ImagSigma = str;
    }
  }
  getline(gout,str); // # number of operators needed
  gout>>Osize;
  gout.ignore(1000,'\n');
  getline(gout,str);
  HF_M.resize(Osize);
  for (int op=0; op<Osize; op++){// Operator for the high frequency moments, called Op in the string above
    if (!gout.good()) cerr<<"Something wrong in reading cix file high frequency expansion"<<endl;
    HF_M[op].resize(N_unique_fl,nsize+1);
    for (int i=1; i<=nsize; i++){
      for (int ib=0; ib<N_unique_fl; ib++){
	int it, size1, size2;
	gout>>it>>size1>>size2;
	if (it!=i) cerr<<"Something wrong reading cix file"<<endl;
	if (size1!=msize_[i] || size2!=msize_[i]) cerr<<"Sizes of HFM are wrong!"<<endl;
	HF_M[op](ib,i).resize(size2,size1); // constructor changes i->ii
	for (int i1=0; i1<size1; i1++){
	  for (int i2=0; i2<size2; i2++){
	    double m;
	    gout>>m;
	    HF_M[op](ib,i)(i2,i1)=m;    
	  }
	}
	gout.ignore(1000,'\n');
      }
    }
    getline(gout,str);
  }
}

void ClusterData::EvaluateMa(double beta, double mu, double U)
{
  Ma.resize(N_flavors);
  for (size_t i=0; i<Ma.size(); i++) Ma[i].resize(size()+1,max_size);
  
  function2D<double> D(max_size,max_size);
  for (int i=1; i<=size(); i++){
    for (int op=0; op<N_flavors; op++){
      int inew = F_i(2*op+1,i);
      if (inew==0) continue;
      D.MProduct(F_M(2*op,inew),F_M(2*op+1,i));
      for (int m=0; m<D.size_N(); m++)	Ma[op](i,m) = D(m,m);
    }
  }

  for (int i=1; i<=nsize; i++){
    double dE = -Ns[i]*mu + 0.5*Ns[i]*(Ns[i]-1)*U;
    for (int m=0; m<Ene[i].size(); m++) Ene(i,m) += dE;
  }

  //  cout<<"GS are:"<<endl;
  gs_ind.resize(nsize+1);
  for (int i=1; i<=nsize; i++){
    int imin = 0;
    for (int m=1; m<Ene[i].size(); m++)
      if (Ene(i,m)<Ene(i,imin)) imin = m;
    gs_ind[i] = imin;
    //    cout<<setw(3)<<i<<" "<<setw(4)<<gs_ind[i]<<endl;
  }
 
  Zatom=0;
  for (int i=1; i<=size(); i++)
    for (int m=0; m<msize(i); m++)
      Zatom += Number(1,-Ene[i][m]*beta);
  
  Patom.resize(size());
  Patom_.resize(size()+1,max_size);
  for (int i=1; i<=size(); i++){
    double sum=0;
    for (int m=0; m<msize(i); m++){
      Patom_(i,m) = divide(Number(1,-Ene[i][m]*beta),Zatom);
      sum += Patom_(i,m);
    }
    Patom[i-1] = sum;
  }

  natom.resize(N_flavors);
  for (int b=0; b<N_flavors; b++){
    double sum=0;
    for (int i=1; i<=size(); i++)
      for (int m=0; m<msize(i); m++)
	sum += Patom_(i,m)*Ma[b](i,m);
    natom[b] = sum;
  }

}

bool CheckStream(istream& inputf, int& n, int& m)
{
  istream input(inputf.rdbuf());
  input.seekg(0,ios::beg);

  n=0;
  string str;
  bool begincomment = false;
  getline(input,str);
  if (!input.good()) {
    cerr << "ERROR: Wrong file format for input stream or no data!" << endl;
    return false;
  }
  if (str.find('#')<string::npos){
    begincomment = true;
  }else n++;
  
  getline(input,str); n++;
  stringstream oneline(str);
  m=0; double t;
  while (oneline){oneline>>t; m++;}
  m--;
  while (input){ getline(input,str); n++;}
  n--;

//   clog << " Number of entries: "<< n <<endl;
//   clog << " Number of columns: "<< m <<endl;

  inputf.seekg(0,ios::beg);
  if (begincomment) getline(inputf, str);
  if (!inputf){ cerr<<"Reopening didn't suceeded!"<<endl; return false;}
  return true;
}

bool ReadData(istream& input, mesh1D& om, function2D<dcomplex>& fi, int mdata, int baths)
{
  clog<<"Reading data with "<<om.size()<<" entries"<<endl;
  if (mdata<2*baths+1) {cerr<<"Not enough columns in input Delta file!"<<endl; return false;}
  vector<double> data(mdata);
  int i=-1;
  while (input && ++i<om.size()){
    for (int j=0; j<mdata; j++) input>>data[j];
    input.ignore(500,'\n');
    for (int j=0; j<baths; j++) fi[j][i] = dcomplex(data[2*j+1],data[2*j+2]);
    om[i] = data[0];
  }
  input.clear();              // forget we hit the end of file
  input.seekg(0, ios::beg);   // move to the start of the file
  string str;
  getline(input, str);        // read comment
  return true;
}

vector<double> FindHighFrequencyExp(int Nf, const mesh1D& omi, const function2D<dcomplex>& fi)
{
  vector<double> ah(fi.size_N());
  for (int b=0; b<fi.size_N(); b++){
    double S=0, Sx=0, Sy=0, Sxx=0, Sxy=0;
    for (int j=omi.size()-Nf; j<omi.size(); j++){
      double y = fi[b][j].imag()*omi[j];
      double x = omi[j];
      Sy += y;
      Sx += 1/(x*x);
      Sxx += 1/(x*x*x*x);
      Sxy += y/(x*x);
      S++;
    }
    double dd = S*Sxx-Sx*Sx;
    double a = (Sxx*Sy-Sx*Sxy)/dd;
    //    double bx = (S*Sxy-Sx*Sy)/dd;
    //    clog<<a<<" "<<bx<<endl;
    ah[b] = -a;
  }
  return ah;
}

void CreateLogMesh(int M1, int M2, double T, const mesh1D& omi, mesh1D& oms)
{
  oms.resize(M1+M2+1);
  
  for (int i=0; i<M1; i++) oms[i] = omi[i];

  double alpha = log((omi.size()-1.)/M1)/(M2-1.);
  
  int inp = static_cast<int>(0.5*(oms[M1-1]/(M_PI*T)-1));
  int l=0;
  for (int i=0; i<M2; i++){
    int in = static_cast<int>(M1*exp(alpha*i)+0.5);
    if (in!=inp) oms[M1+(l++)] = (2*in+1)*M_PI*T;
    inp = in;
  }
  int last = static_cast<int>(0.5*(omi.last()/(M_PI*T)-1));
  if (inp!=last) oms[M1+(l++)] = (2*last+1)*M_PI*T;
  oms.resize(l);
  oms.SetUp(0.0);
}


void InverseFourier(const function<dcomplex>& Giom, const mesh& iom, function<double>& Gtau, const mesh& tau, double ah, double beta)
{
  for (int t=0; t<tau.size(); t++){
    double sum=0;
    for (int n=0; n<iom.size(); n++)
      sum += cos(iom[n]*tau[t])*Giom[n].real() + sin(iom[n]*tau[t])*(Giom[n].imag()+ah/iom[n]);
    Gtau[t] = 2*sum/beta-0.5*ah;
  }
  // Correction 1/omega^2 (sometimes usefull)
  // corrected 16.12.2006
  double df = Giom.last().real()*iom.last()/M_PI;
  Gtau[0] += df;
  Gtau[Gtau.size()-1] -= df;  
}

bool ReadDelta(int Ns, istream& input, int n, int m, mesh1D& omi, function2D<dcomplex>& Deltaw, vector<double>& ah, double beta)
{
  mesh1D tomi(n);
  function2D<dcomplex> tDeltaw(Ns,tomi.size());
  if (!ReadData(input, tomi, tDeltaw, m, Ns)) return false;
  tomi.SetUp(0);
  ah = FindHighFrequencyExp(500, tomi, tDeltaw);

  omi.resize(tomi.size());
  Deltaw.resize(Ns,omi.size());
  for (int ib=0; ib<Ns; ib++)
    for (int i=0; i<omi.size(); i++){
      omi[i] = (2*i+1)*M_PI/beta;
      if (omi[i]<tomi.last()+M_PI/beta) // correction 16.12.2006
	Deltaw[ib][i] = tDeltaw[ib](tomi.Interp(omi[i]));
      else
	Deltaw[ib][i] = dcomplex(0,-ah[ib]/omi[i]);
    }

  // This is exact for half-filling non-interacting case
//   for (int i=0; i<omi.size(); i++)
//     Deltaw[0][i] = dcomplex(0,0.5*omi[i]*(1-sqrt(1+1/(omi[i]*omi[i]))));

  return true;
}

void DeltaFourier(int Ntau, double beta, mesh1D& tau, vector<spline1D<double> >& Delta, const mesh1D& omi, const function2D<dcomplex>& Deltaw, const vector<double>& ah)
{
  tau.MakeEquidistantMesh(Ntau,0,beta);
  for (size_t ib=0; ib<Delta.size(); ib++) Delta[ib].resize(tau.size());
  for (size_t ib=0; ib<Delta.size(); ib++){
    InverseFourier(Deltaw[ib], omi, Delta[ib], tau, ah[ib], beta); // Here needs to be corrected!

    ofstream check(NameOfFile("check.dat",ib).c_str());
    check<<"# ah = "<<ah[ib]<<endl; 
    for (int i=0; i<tau.size(); i++) check<<tau[i]<<" "<<Delta[ib][i]<<endl;

    
    double df0 = (Delta[ib][1]-Delta[ib][0])/(tau[1]-tau[0]);      // If Delta[ib][itau]>0 set it to something negative
    int n = tau.size();
    double dfn = (Delta[ib][n-1]-Delta[ib][n-2])/(tau[n-1]-tau[n-2]);
    Delta[ib].splineIt(tau,df0,dfn);
  }
}

void ClusterData::RecomputeCix(function2D<double>& AProb, int Naver, double treshold, function1D<NState>& praStates)
{

//    ifstream inp("Probability.dat");
//    for (int i=0; i<nsize; i++){
//      for (int j=0; j<msize_[i+1]; j++){
//        int it, mt;
//        double Pi;
//        inp>>it>>mt>>Pi;
//        inp.ignore(500,'\n');
//        if (it-1!=i || mt!=j) cerr<<"Something wromg in reading Probability.dat"<<endl;
//        AProb[i][j] = Pi;
//      }
//    }
//    Naver=1;
  
  for (int i=0; i<nsize; i++) for (int j=0; j<msize_[i+1]; j++) AProb[i][j]/=Naver;
  
  ofstream out("Probability.dat"); out.precision(16);
  out<<"# Naver="<<Naver<<endl;
  for (int j=0; j<AProb.size_N(); j++)
    for (int k=0; k<AProb[j].size(); k++)
      out<<setw(3)<<j+1<<" "<<setw(3)<<k<<" "<<setw(20)<<AProb[j][k]<<endl;
  
  vector<deque<int> > redundant(nsize);
  function1D<int> new_msize_(msize_);
  //  deque<int> remove;
  for (int ist=0; ist<nsize; ist++){
    int ii = ist+1;
    deque<int> redun;
    for (int m=0; m<msize_[ii]; m++)
      if (fabs(AProb(ist,m))<treshold) redun.push_back(m); // this state should be eliminated
    redundant[ist] = redun;
    //    if (redun.size()==msize_[ist+1])  remove.push_back(ist);
    new_msize_[ii] = msize_[ii] - redun.size();
  }

  for (int ist=0; ist<nsize; ist++){
    int ii = ist+1;
    if (new_msize_[ii]==0) praStates[ist].SetIstate(0);
    
    for (int ib=0; ib<N_flavors; ib++){
      int jj = F_i(2*ib,ii);
      if (jj==0) continue;
      int size1 = msize_[ii];
      int size2 = msize_[jj];
      
      for (int k=redundant[ii-1].size()-1; k>=0; k--){
	int ik = redundant[ii-1][k];
	for (int m1=0; m1<size2; m1++)
	  for (int l=ik+1; l<size1; l++){
	    F_M(2*ib,ii)(m1,l-1) = F_M(2*ib,ii)(m1,l);
	    F_M(2*ib+1,jj)(l-1,m1) = F_M(2*ib+1,jj)(l,m1);
	  }
	size1--;
      }
      for (int k=redundant[jj-1].size()-1; k>=0; k--){
	int ik = redundant[jj-1][k];
	for (int m2=0; m2<size1; m2++)
	  for (int l=ik+1; l<size2; l++){
	    F_M(2*ib,ii)(l-1,m2) = F_M(2*ib,ii)(l,m2);
	    F_M(2*ib+1,jj)(m2,l-1) = F_M(2*ib+1,jj)(m2,l);
	  }
	size2--;
      }
      if (size1==0 || size2==0){
	F_i[2*ib][ii]=0;
	F_i[2*ib+1][jj]=0;
	size1=0;
	size2=0;
      }
      double sum=0;
      for (int m1=0; m1<size2; m1++)
	for (int m2=0; m2<size1; m2++)
	  sum += fabs(F_M(2*ib,ii)(m1,m2));
      if (sum<treshold){
	F_i[2*ib][ii]=0;
	F_i[2*ib+1][jj]=0;	
	size1=0;
	size2=0;
      }
      F_M(2*ib,ii).resize(size2,size1);
      F_M(2*ib+1,jj).resize(size1,size2);
    }
  }
  
  for (int ist=0; ist<nsize; ist++){
    int ii = ist+1;
    for (int k=redundant[ii-1].size()-1; k>=0; k--){
      int ik = redundant[ii-1][k];
      for (int l=ik+1; l<msize_[ii]; l++){
	Ene0[ii][l-1] = Ene0[ii][l];
	Spin[ii][l-1] = Spin[ii][l];
      }	
    }
    Ene0[ii].resize(new_msize_[ii]);
    Spin[ii].resize(new_msize_[ii]);
  }
  
  for (int op=0; op<Osize; op++){
    for (int ist=0; ist<nsize; ist++){
      int ii = ist+1;
      for (int ib=0; ib<N_unique_fl; ib++){
	function2D<double>& hf = HF_M[op](ib,ii);
	int size = msize_[ii];
	for (int k=redundant[ii-1].size()-1; k>=0; k--){
	  int ik = redundant[ii-1][k];
	  for (int m1=0; m1<size; m1++) for (int l=ik+1; l<size; l++) hf(m1,l-1) = hf(m1,l);
	  for (int m1=0; m1<size; m1++) for (int l=ik+1; l<size; l++) hf(l-1,m1) = hf(l,m1);
	}
	hf.resize(new_msize_[ii],new_msize_[ii]);
      }
    }
  }
  msize_ = new_msize_;
  
  ofstream gout("new.cix"); gout.precision(12);
  gout<<"# Updated Cix file for cluster DMFT with CTQMC"<<endl;
  gout<<"# cluster_size, number of states, number of baths, maximum_matrix_size"<<endl;
  gout<<Nm<<" "<<nsize<<" "<<N_ifl<<" "<<max_size<<endl;
  gout<<"# baths, dimension, symmetry"<<endl;
  for (int i=0; i<N_ifl; i++){
    gout<<left<<setw(3)<<i<<" "<<setw(4)<<ifl_dim[i]<<" ";
    for (int b=0; b<sqr(ifl_dim[i]); b++){
      if (sign[i][b]<0) gout<<"-";
      gout<<bfl_index[i][b];
      if (conjg[i][b]) gout<<"*";
      gout<<" ";
    }
    gout<<"   "<<setw(3)<<gflip_index[i]<<endl;
  }
  gout<<"# cluster energies for non-equivalent baths, eps[k]"<<endl;

  for (int ik=0; ik<epsk.size(); ik++) gout<<epsk[ik]<<" ";
  gout<<endl;
  
  gout<<right;
  gout<<" #   N   K   Sz size"<<endl;
  for (int ist=0; ist<nsize; ist++){
    int ii = ist+1;
    gout<<setw(2)<<ii<<" "<<setw(3)<<Ns[ii]<<" "<<setw(3)<<Ks[ii]<<" "<<setw(4)<<right<<Sz[ii]<<" "<<setw(3)<<msize_[ii]<<"    ";
    for (int ib=0; ib<N_flavors; ib++) gout<<setw(3)<<F_i(2*ib,ist+1)<<" ";
    //    gout<<" ";
    
    for (int is=0; is<Ene0[ii].size(); is++) gout<<setw(18)<<Ene0(ii,is)<<" ";
    for (int is=0; is<Spin[ii].size(); is++) gout<<setw(18)<<Spin(ii,is)<<" ";
    gout<<endl;
  }
  gout<<"# matrix elements"<<endl;
  for (int ist=0; ist<nsize; ist++){
    for (int ib=0; ib<N_flavors; ib++){
      gout<<setw(2)<<ist+1<<" "<<setw(3)<<F_i(2*ib,ist+1)<<"  ";
      function2D<double>& fm = F_M(2*ib,ist+1);
      gout<<setw(3)<<fm.size_Nd()<<" "<<setw(3)<<fm.size_N()<<"  ";
      for (int im1=0; im1<fm.size_Nd(); im1++)
	for (int im2=0; im2<fm.size_N(); im2++)
	  gout<<setw(20)<<fm[im2][im1]<<" ";
      gout<<endl;
    }
  }
  gout<<"# high frequency expansion (Real and Imaginary part of Sigma)"<<endl;
  gout<<RealSigma<<endl;
  gout<<ImagSigma<<endl;
  gout<<"# number of operators needed"<<endl;
  
  gout<<Osize<<endl;
  for (int op=0; op<Osize; op++){
    gout<<"# high frequency moment matrix elements, called Op"<<op+1<<" above"<<endl;
    for (int i=1; i<=nsize; i++){
      for (int ib=0; ib<N_unique_fl; ib++){
	int size = msize_[i];
	gout<<setw(2)<<i<<"  "<<setw(3)<<size<<" "<<setw(3)<<size<<"      ";
	function2D<double>& hf = HF_M[op](ib,i);
	for (int i1=0; i1<size; i1++)
	  for (int i2=0; i2<size; i2++)
	    gout<<setw(20)<<hf(i2,i1)<<" ";
	gout<<endl;
      }
    }
  }

  Ene = Ene0;
  for (int i=1; i<=nsize; i++){
    double dE = -Ns[i]*common::mu + 0.5*Ns[i]*(Ns[i]-1)*common::U;
    for (int m=0; m<Ene[i].size(); m++) Ene(i,m) += dE;
  }
}

void Inverse_Gf(int dim, int N_unique_fl, const deque<int>& bfl_index, const deque<int>& sign, const deque<int>& conjg, int im, const function2D<dcomplex>& Gf, function2D<dcomplex>& Gf_1)
{
  static function2D<dcomplex> g(dim,dim), gi(dim,dim);
  if (dim==1){
    Gf_1[bfl_index[0]][im] += 1/Gf[bfl_index[0]][im];
    return;
  }
  g.resize(dim,dim);
  gi.resize(dim,dim);
  int ib=0;
  for (int i1=0; i1<dim; i1++){
    for (int i2=0; i2<dim; i2++){
      dcomplex gf = Gf(bfl_index[ib],im);
      if (sign[ib]==-1) gf*=-1;
      if (conjg[ib]) gf = gf.conj();	  
      g[i1][i2] = gf;
      ib++;
    }
  }
  gi.Inverse(g);
  ib=0;
  for (int i1=0; i1<dim; i1++){
    for (int i2=0; i2<dim; i2++){
      dcomplex gf = gi[i1][i2];
      if (sign[ib]==-1) gf*=-1;
      if (conjg[ib]) gf = gf.conj();
      Gf_1(bfl_index[ib], im) += gf;
      ib++;
    }
  }
}