void StartStateEvolution(function2D<NState>& state_evolution_left,
			 function2D<NState>& state_evolution_right,
			 const function1D<NState>& praStates,
			 NOperators& Operators, const ClusterData& cluster,
			 function1D<Number>& Trace, function2D<double>& Prob)
{
  if (!Operators.empty()){cerr<<"Operators are not empty!"<<endl;}
  function2D<Number> Projection(Prob.size_N(),common::max_size);
  NState nstate;
  Number ms=0;
  for (int i=0; i<cluster.size(); i++){
    if (praStates[i].empty()) {
      state_evolution_left[i].resize(0);
      state_evolution_right[i].resize(0);
      continue;
    }
    state_evolution_left[i].resize(1);
    state_evolution_right[i].resize(1);
    nstate = praStates[i];
    state_evolution_left(i,0).Evolve(nstate, cluster,  Operators.exp_(0));
    state_evolution_right(i,0).Evolve(nstate, cluster, Operators.exp_(0));
    Trace[i] = state_evolution_left[i][0].Project_to_Pra(praStates[i], Projection[i]);//TraceProject(praStates[i]);
    ms += Trace[i];
  }
  for (int i=0; i<cluster.size(); i++)
    for (int l=0; l<Prob[i].size(); l++)
      Prob[i][l] = divide(Projection[i][l],ms);
}

Number UpdateStateEvolution(function2D<NState>& state_evolution_left, function2D<NState>& state_evolution_right,
			    int op_i, int op_f, NOperators& Operators, const function1D<NState>& praStates, const ClusterData& cluster,
			    function1D<Number>& Trace, function2D<double>& Prob)
{
  static function2D<Number> Projection(Prob.size_N(),common::max_size);
  static NState nstate, mstate;
  Number ms=0;
  Projection=0;
  for (int ist=0; ist<praStates.size(); ist++){// op_i is the earliest operator inserted - from this point on, the evolution needs to be changed
    Trace[ist]=0;
    if (op_i>0 && (state_evolution_left[ist].size()<op_i || state_evolution_left[ist][op_i-1].empty())) continue;// this product matrix element is zero
    int ip = op_i;
    if (ip==0) nstate = praStates[ist];// we need to start from scratch because op_i iz zero
    else nstate = state_evolution_left[ist][ip-1];// we start from the previous operator and update from there on
    if (nstate.empty()) continue;
      
    bool empty=false;
    
    for (; ip<Operators.size(); ip++){        // go over the rest of the operators
      mstate.Evolve(nstate, cluster, Operators.exp_(ip)); // state evolution due to local Hamiltonian (is diagonal)
      int iop = Operators.typ(ip);
      nstate.apply(cluster.FM(iop), cluster.Fi(iop), mstate);        // the next operator is applied
      if (nstate.empty()) {empty=true; break;}// maybe matrix element is zero
      state_evolution_left[ist][ip] = nstate;      // it is not, store it for next time
    }

    if (!empty){ // from the last operator to beta
      mstate.Evolve(nstate, cluster, Operators.exp_last()); // The last evlution of the state due to H_loc.
      state_evolution_left[ist][ip++] = mstate; // store it
    }
    if (ip>state_evolution_left.fullsize_Nd()) {cerr<<"Trying to resize to "<<ip<<" op_size="<<Operators.size()<<endl;}
    state_evolution_left[ist].resize(ip); // we have stored state evolution, remember the new size
    Number mm=0;
    if (!empty) mm = mstate.Project_to_Pra(praStates[ist], Projection[ist]);
    Trace[ist] = mm;
    ms += mm;
  }

  Number ms_right=0;
  for (int ist=0; ist<praStates.size(); ist++){// op_i is the earliest operator inserted - from this point on, the evolution needs to be changed
    int nsize = Operators.size();
    int ip = nsize-1-op_f;
    if (ip>0 && (state_evolution_right[ist].size()<ip || state_evolution_right[ist][ip-1].empty())) continue;// this product matrix element is zero
    if (ip==0) nstate = praStates[ist];// we need to start from scratch because op_i iz zero
    else nstate = state_evolution_right[ist][ip-1];// we start from the previous operator and update from there on
    if (nstate.empty()) continue;
    
    bool empty=false;

    if (!empty && ip<nsize){ // from the last operator to beta
      if (ip>0) mstate.Evolve(nstate, cluster, Operators.exp_(nsize-ip)); // The last evlution of the state due to H_loc.
      else 	mstate.Evolve(nstate, cluster, Operators.exp_last()); // The last evlution of the state due to H_loc.
      int iop = Operators.typ_transpose(nsize-1-ip);// HERE iop==-1 !!!!
      nstate.apply(cluster.FM(iop), cluster.Fi(iop), mstate);        // the next operator is applied
      if (nstate.empty()) {empty=true;}// maybe matrix element is zero
      if (!empty) state_evolution_right[ist][ip++] = nstate; // store it
    }
    if (!empty){
      for (; ip<Operators.size(); ip++){        // go over the rest of the operators
	mstate.Evolve(nstate, cluster, Operators.exp_(nsize-ip)); // state evolution due to local Hamiltonian (is diagonal)
	int iop = Operators.typ_transpose(nsize-1-ip);
	nstate.apply(cluster.FM(iop), cluster.Fi(iop), mstate);        // the next operator is applied
	if (nstate.empty()) {empty=true; break;}// maybe matrix element is zero
	state_evolution_right[ist][ip] = nstate;      // it is not, store it for next time
      }
    }
    if (!empty){
      mstate.Evolve(nstate, cluster, Operators.exp_(0)); // The last evlution of the state due to H_loc.
      state_evolution_right[ist][ip++] = mstate; // store it
    }
    if (ip>state_evolution_right.fullsize_Nd()) {cerr<<"Trying to resize to "<<ip<<" op_size="<<Operators.size()<<endl;}
    state_evolution_right[ist].resize(ip); // we have stored state evolution, remember the new size
    Number mm=0;
    if (!empty) {
      mm = mstate.TraceProject(praStates[ist]); // need the trace. The last state needs to be projected to the first state
    }
    ms_right += mm;
  }

  if (fabs(fabs(divide(ms_right,ms))-1)>1e-6){
    cerr<<"ms left and ms right are not the same "<<ms<<" "<<ms_right<<endl;
    if (fabs(divide(ms_right,ms))>1) ms = ms_right;
  }
  double sum=0;
  for (int ist=0; ist<praStates.size(); ist++){
    for (int l=0; l<Prob[ist].size(); l++){
      Prob[ist][l] = divide(Projection[ist][l],ms);
      sum += Prob[ist][l];
    }
  }

  if (common::PreciseP){ // computes probability more precisely
    int nsize = Operators.size();
    if (nsize==0){
      for (int ist=0; ist<praStates.size(); ist++)
	for (int m=0; m<Prob[ist].size(); m++)
	  Prob(ist,m) = cluster.P_atom(ist,m);
    }else{
      Prob=0;
      for (int ist=0; ist<praStates.size(); ist++){
	if (state_evolution_left[ist].size()<=Operators.size() || state_evolution_right[ist].size()<=Operators.size()) continue;
	if (op_i>0 && (state_evolution_left[ist].size()<op_i || state_evolution_left[ist][op_i-1].empty())) continue;// this product matrix element is zero
      
	for (int l=0; l<nsize-1; l++){
	  NState& lstate = state_evolution_left(ist,l);
	  NState& rstate = state_evolution_right(ist,nsize-l-2);
	
	  const function<double>& expn = Operators.exp_(l+1)[lstate.istate];
	  function<double>& Prb = Prob[lstate.istate-1];
	  double dtau = (Operators.t(l+1)-Operators.t(l))/common::beta;
	  double dtau_ms_mantisa = dtau/ms.mantisa;
	  double out_exponents = lstate.exponent+rstate.exponent-ms.exponent;
	  
	  for (int n=0; n<lstate.M.size_N(); n++){
	    double mantisa=0;
	    for (int m=0; m<lstate.M.size_Nd(); m++) mantisa += lstate.M(n,m)*rstate.M(n,m);
	    Prb[n] += mantisa*dtau_ms_mantisa*exp(expn[n]+out_exponents);
	  }
	}
    
	NState& lstate = state_evolution_left(ist,nsize);
	function<double>& Prb = Prob[lstate.istate-1];
    
	double dtau = 1 - (Operators.t(nsize-1)-Operators.t(0))/common::beta;
	for (int n=0; n<lstate.M.size_N(); n++)
	  Prb[n] += lstate.M(n,n)/ms.mantisa*exp(lstate.exponent-ms.exponent)*dtau;
      }
    }
    double sum=0;
    for (int ist=0; ist<praStates.size(); ist++)
      for (int l=0; l<Prob[ist].size(); l++)
	sum += Prob[ist][l];
    
    //    cout<<nsize<<" "<<sum<<endl;
  }

  return ms*Operators.sign();

}

Number ComputeTrace(const function2D<NState>& state_evolution_left,
		    const function2D<NState>& state_evolution_right,
		    int op_i, int op_f,
		    const NOperators& Operators, const function1D<NState>& praStates,
		    const ClusterData& cluster,
		    int dsize)
{
  
  static NState nstate, mstate;
  Number ms_new=0;
  for (int ist=0; ist<praStates.size(); ist++){
    if (op_i>0 && (state_evolution_left[ist].size()<op_i || state_evolution_left[ist][op_i-1].empty())) continue;
    int ip = op_i;
    if (ip==0) nstate = praStates[ist];
    else nstate = state_evolution_left[ist][ip-1];
    if (nstate.empty()) continue;
      
    Number mm_new=0;
    bool empty=false;
    int nsize = Operators.size();
    
    for (; ip<=op_f; ip++){
      mstate.Evolve(nstate, cluster, Operators.tr_exp_(ip));
      int iop = Operators.tr_typ(ip);
      nstate.apply(cluster.FM(iop), cluster.Fi(iop), mstate);
      if (nstate.empty()) {empty=true; break;}
    }
    if (!empty){
      if (op_f<nsize+dsize-1){
	if (state_evolution_right[ist].size()>nsize+dsize-op_f-2 &&
	    nstate.istate==state_evolution_right[ist][nsize+dsize-op_f-2].istate){
	  mstate.Evolve(nstate, cluster, Operators.tr_exp_(ip));
	  mm_new = mstate.ScalarProduct(state_evolution_right[ist][nsize+dsize-op_f-2]);
	}else{
	  mm_new = 0;
	}
      }else{
	mstate.Evolve(nstate, cluster, Operators.tr_exp_last());
	mm_new = mstate.TraceProject(praStates[ist]);
      }
    }
    ms_new += mm_new;
  }
  return abs(ms_new);
}

Number ComputeGlobalTrace(const NOperators& Operators, const function1D<NState>& praStates, const ClusterData& cluster, const function<int>& gflip_fl)
{
  static NState nstate, mstate;
  Number ms_new=0;
  int nsize = Operators.size();
  for (int ist=0; ist<praStates.size(); ist++){
    int ip = 0;
    nstate = praStates[ist];
    if (nstate.empty()) continue;
    
    Number mm_new=0;
    bool empty=false;
    
    for (; ip<nsize; ip++){
      mstate.Evolve(nstate, cluster, Operators.exp_(ip));
      int iop = Operators.typ(ip);
      int nfl = gflip_fl[iop/2];
      int nop = 2*nfl + iop%2;
      
      nstate.apply(cluster.FM(nop), cluster.Fi(nop), mstate);
      if (nstate.empty()) {empty=true; break;}
    }
    if (!empty){
      mstate.Evolve(nstate, cluster, Operators.exp_last());
      mm_new = mstate.TraceProject(praStates[ist]);
    }
    ms_new += mm_new;
  }
  return abs(ms_new);
}

void ComputeMoments(const ClusterData& cluster, const mesh1D& omi, int iom_start, double nf,
		    const function1D<double>& mom, const function1D<double>& kaver, const function2D<double>& nt,
		    function2D<dcomplex>& Sigma, bool CorrectMoments=true, int Ncorrect=-1, int aom=3, double sderiv=0.1)
{
  if (Ncorrect<0) Ncorrect=cluster.N_unique_fl;
  stringstream script; script.precision(15);
  script<<"#!/usr/bin/perl\n";
  script<<"use Math::Trig;\nuse Math::Complex;\n\n"; 
  script<<"$nf="<<nf<<";\n";
  script<<"$T="<<(1/common::beta)<<";\n";
  script<<"$U="<<common::U<<";\n";
  script<<"$mu="<<common::mu<<";\n";
  script<<"$Nm="<<cluster.Nm<<";\n";
  script<<"$NK="<<cluster.N_unique_fl<<";\n";
  script<<"$NKc="<<Ncorrect<<";\n";
  for (int i=0; i<nt.size_N(); i++){
    script<<"@nt"<<i<<"=(";
    for (int j=0; j<nt.size_Nd()-1; j++) script<<nt[i][j]<<",";
    script<<nt[i].last()<<");\n";
  }
  for (int op=0; op<cluster.HF_M.size(); op++){
    script<<"@Op"<<op+1<<"=(";
    for (int ik=0; ik<cluster.N_unique_fl; ik++){
      script<<mom[ik+op*cluster.N_unique_fl];
      if (ik<cluster.N_unique_fl-1) script<<",";
      else script<<");\n";
    }
  }
  script<<"@epsk=(";
  for (int ik=0; ik<cluster.N_unique_fl; ik++){
    script<<cluster.epsk[ik];
    if (ik<cluster.N_unique_fl-1) script<<",";
    else script<<");\n";
  }
  
  script<<"@kaver=(";
  for (int ik=0; ik<kaver.size(); ik++){
    script<<kaver[ik];
    if (ik<kaver.size()-1) script<<",";
    else script<<");\n";
  }
  
  script<<"@om=(";
  for (int im=iom_start; im<omi.size(); im++){
    script<<omi[im];
    if (im<omi.size()-1) script<<",";
    else script<<");\n";
  }
  script<<"\n";
  script<<"for ($im=0; $im<=$#om; $im++){\n";
  script<<"  for ($k=0; $k<$NK; $k++){\n";
  script<<"    "<<cluster.RealSigma<<";"<<endl;
  script<<"    "<<cluster.ImagSigma<<";"<<endl;
  script<<"    $$Sig[$k][$im] = $ReSigma + $ImSigma*i;\n";
  script<<"  }\n";
  //  script<<"  print \"\\n\";\n";
  script<<"}\n";

  if (CorrectMoments){
    // Corrects high frequency by interpolation
    script<<endl<<endl;
    script<<"@sb=(";
    for (int ik=0; ik<cluster.N_unique_fl; ik++){
      dcomplex sb=0;
      for (int im=iom_start-aom; im<iom_start; im++) sb += Sigma[ik][im];
      sb/=aom;
      script<<sb.real()<<" ";
      if (sb.imag()>=0) script<<"+";
      script<<sb.imag()<<"*i";
      if (ik<cluster.N_unique_fl-1) script<<",";
      else script<<");"<<endl;
    }
    script<<endl;
    script<<"$omb = "<<omi[iom_start-aom/2]<<";"<<endl;
    script<<"$sderiv = "<<sderiv<<";"<<endl<<endl;
    script<<"for ($k=0; $k<$NKc; $k++) {$se[$k] = Re($$Sig[$k][$#om-1]);}"<<endl;
  
    script<<"for ($im=0; $im<=$#om; $im++){"<<endl;
    script<<"    for ($k=0; $k<$NKc; $k++){"<<endl;
    script<<"	$sr = $se[$k] + (Re($sb[$k])-$se[$k])*($omb/$om[$im])**2;"<<endl;
    script<<"	$$Sig[$k][$im] = $sr + Im($$Sig[$k][$im])*i;"<<endl;
    script<<"    }"<<endl;
    script<<"}"<<endl;
    script<<endl;
    script<<"for ($k=0; $k<$NKc; $k++){"<<endl;
    script<<"    for ($iw=1; $iw<$#om; $iw++){"<<endl;
    script<<"	$df0 = Im($$Sig[$k][$iw]-$sb[$k])/($om[$iw]-$omb);"<<endl;
    script<<"	$df1 = Im($$Sig[$k][$iw]-$$Sig[$k][$iw-1])/($om[$iw]-$om[$iw-1]);"<<endl;
    script<<"	if (abs($df0-$df1)<$sderiv) {last;}"<<endl;
    script<<"    }"<<endl;
    script<<"    $niw[$k] = $iw;"<<endl;
    script<<"    $se[$k] = Im($$Sig[$k][$iw]);"<<endl;
    script<<"    $oe[$k] = $om[$iw];"<<endl;
    script<<"}"<<endl;
    script<<endl;
    script<<"for ($k=0; $k<$NKc; $k++){"<<endl;
    script<<"    for ($im=0; $im<$niw[$k]; $im++){"<<endl;
    script<<"	$si = Im($sb[$k]) + ($se[$k]-Im($sb[$k]))*($om[$im]-$omb)/($oe[$k]-$omb);"<<endl;
    script<<"	$$Sig[$k][$im] = Re($$Sig[$k][$im]) + $si*i;"<<endl;
    script<<"    }"<<endl;
    script<<"}"<<endl;
  }
  script<<endl;
  script<<"for ($im=0; $im<=$#om; $im++){"<<endl;
  script<<"  print $om[$im], \"  \";"<<endl;
  script<<"  for ($k=0; $k<$NK; $k++){"<<endl;
  script<<"    print Re($$Sig[$k][$im]), \" \", Im($$Sig[$k][$im]), \"  \";"<<endl;
  script<<"  }\n";
  script<<"  print \"\\n\";\n";
  script<<"}\n";
  script<<endl;
  
  
  ofstream scrpt("pscript.pl");
  scrpt<<script.str()<<endl;
  
  string command;
  command = "perl pscript.pl > moments.temp";
  int ifail = system(command.c_str());

  ifstream inp("moments.temp");
  if (ifail){cerr<<"Seems your moments are not working or could not find perl!"<<endl;}
  for (int im=iom_start; im<omi.size(); im++){
    double omega;
    inp>>omega;
    if (fabs(omega-omi[im])>1e-6) {cerr<<"Did not succeed in reading moments. Exiting!"<<endl; break;}
    for (int ik=0; ik<cluster.N_unique_fl; ik++) inp>>Sigma(ik,im);
    inp.ignore(1000,'\n');
    if (!inp) {cerr<<"Did not succeed in reading moments. Exiting!"<<endl; break;}
  }
}