C++ lst::begin方法代码示例

UFXmap UF(lst k_lst,lst p_lst,lst subs_lst, unsigned int displacement )
{
  /// constructing expr Xi*Pi
  ex fprop = 0;
  ex sDtmp;
  string str;
  lst xp;
  for(lst::const_iterator it = p_lst.begin();it!=p_lst.end();++it)
    {
      str = "x_" + boost::lexical_cast<string>(displacement+std::distance(p_lst.begin(),it));
      xp.append(get_symbol(str)); // storing list of X identities
      fprop += get_symbol(str) * (*it); 
    }
  fprop = fprop.expand();
  sDtmp = fprop;
  //cout<<fprop<<endl;
  matrix M(k_lst.nops(),k_lst.nops());
  matrix Q(k_lst.nops(),1);//column
  GiNaC::numeric half(1,2);
  for(lst::const_iterator itr = k_lst.begin();itr!=k_lst.end();++itr)
    for(lst::const_iterator itc = itr;itc!=k_lst.end();++itc)
      if(itr == itc)
        {
          M(distance(k_lst.begin(),itr),distance(k_lst.begin(),itc)) = -1*fprop.coeff((*itr)*(*itc),1);
          //cout<<(*itr)*(*itc)<<"M("<<distance(k_lst.begin(),itr)<<","<<distance( k_lst.begin(),itc)<<") coeff "<<fprop.coeff((*itr)*(*itc),1)<<endl;
          sDtmp -= (*itr)*(*itc)*fprop.coeff((*itr)*(*itc),1);
        }
      else
        {
          M(distance(k_lst.begin(),itr),distance( k_lst.begin(),itc)) = -1*half*fprop.coeff((*itr),1).coeff((*itc),1);
          M(distance(k_lst.begin(),itc),distance(k_lst.begin(),itr)) = -1*half*fprop.coeff((*itr),1).coeff((*itc),1);
          //cout<<(*itr)*(*itc)<<"M("<<distance( k_lst.begin(),itr)<<","<<distance( k_lst.begin(),itc)<<") coeff "<<fprop.coeff((*itr),1).coeff((*itc),1)<<endl;
          sDtmp -= (*itr)*(*itc)*fprop.coeff((*itr),1).coeff((*itc),1);
        }
  cout<<"M: "<<M<<endl;
  sDtmp = sDtmp.expand();
  //cout<<"Expr linear on external momentum: "<<sDtmp.expand()<<endl;

  for(lst::const_iterator itr = k_lst.begin();itr!=k_lst.end();++itr)
    {
      Q(distance( k_lst.begin(),itr),0) = half*sDtmp.coeff((*itr),1);
      sDtmp -= (*itr)*sDtmp.coeff((*itr),1);
    }
  //  cout<<"Q: "<<Q<<endl;
  sDtmp = sDtmp.expand();
  ex minusJ = sDtmp;
   cout<<"-J: "<<minusJ<<endl;
  ex U = M.determinant();
  ex F = expand(M.determinant()*(minusJ+Q.transpose().mul(M.inverse()).mul(Q)(0,0)));
  lst lp;
  F=F.normal();
  cout<<"U= "<<U<<endl<<"F= "<<F<<endl;
  //  cout<<"pol_list "<<lp<<endl;
  U=U.subs(subs_lst,subs_options::algebraic);
  F=F.subs(subs_lst,subs_options::algebraic);
  cout<<"ALGEBRAIC: U= "<<U<<endl<<"F= "<<F<<endl;
  //  cout<<"UF_WORK"<<endl;
  return fusion::make_map<UFX::U,UFX::F,UFX::xlst>(U,F,xp);
}

bool zero_volume(const lst& pole_list,const lst& w_list)
{
  try{
  using namespace lemon;
  GlpkLp lp;
  exhashmap<LpBase::Col> col_map;
  for(lst::const_iterator wi = w_list.begin();wi!=w_list.end();++wi)
    {
      col_map[*wi] = lp.addCol();
    }
  for(lst::const_iterator pit = pole_list.begin();pit!= pole_list.end();++pit)
    {
      ex tmp_expr = *pit;
      Lp::Expr constr_expr;
      for(lst::const_iterator wi = w_list.begin();wi!=w_list.end();++wi)
	{
	  
	  ex wi_coeff = tmp_expr.coeff(*wi);
	  tmp_expr-=(*wi)*wi_coeff;
	  if(is_a<numeric>(wi_coeff))
	    constr_expr+=ex_to<numeric>(wi_coeff).to_double()*col_map[*wi];
	  else throw std::logic_error(std::string("Non numeric coefficient in pole term. "));
	}
      if(is_a<numeric>(tmp_expr))
	lp.addRow(-ex_to<numeric>(tmp_expr).to_double(),constr_expr,Lp::INF);
      else 
	{
	  cout<< tmp_expr<<endl;
	  throw std::logic_error(std::string("Lower bound is not a numeric"));
	}
    }
    
  double l_bound,u_bound;
  for(lst::const_iterator wi = w_list.begin();wi!=w_list.end();++wi)
    {
      lp.min();
      lp.obj(col_map[*wi]);
      lp.solve();
      if (lp.primalType() == Lp::OPTIMAL) 
	l_bound = lp.primal();
      else return true;
      lp.max();
      lp.obj(col_map[*wi]);
      lp.solve();
      if (lp.primalType() == Lp::OPTIMAL) 
	u_bound = lp.primal();
      else return true;
      cout<<"ub: "<<u_bound<<" lb: "<<l_bound<<endl;
      if((u_bound - l_bound) <= 0)
        return true;
    }
  return false;
  }catch(std::exception &p)
    {
      throw std::logic_error(std::string("In function \"zero_volume\":\n |___> ")+p.what());
    }
}

//
// Replace each var in e with delay(var,lag).
//
ex delay_vars(const ex& e, const ex& lag, const lst& vars)
{
    ex g = e;
    for (lst::const_iterator i = vars.begin(); i != vars.end(); ++i) {
        g = g.subs(*i == delay(*i, lag));
    }
    return g;
}

exset lst2set(const lst & l) {
  exset s;
  lst::const_iterator it = l.begin();
  while (it != l.end()) {
    s.insert(*it);
    ++it;
  }
  return s;
}

/** Substitute objects in an expression (syntactic substitution) and return
 *  the result as a new expression. */
ex ex::subs(const lst & ls, const lst & lr, unsigned options) const
{
	GINAC_ASSERT(ls.nops() == lr.nops());

	// Convert the lists to a map
	exmap m;
	for (lst::const_iterator its = ls.begin(), itr = lr.begin(); its != ls.end(); ++its, ++itr) {
		m.insert(std::make_pair(*its, *itr));

		// Search for products and powers in the expressions to be substituted
		// (for an optimization in expairseq::subs())
		if (is_exactly_a<mul>(*its) || is_exactly_a<power>(*its))
			options |= subs_options::pattern_is_product;
	}
	if (!(options & subs_options::pattern_is_product))
		options |= subs_options::pattern_is_not_product;

	return bp->subs(m, options);
}

//
// Merge the expressions in exprs and formulas into a single expression.
// They are "merged" using the equality relation.  That is, if
// exprs = {e1, e2, e3} and formulas = {f1, f2}, the result is
// e1 == (e2 == (e3 == (f1 == f2)))
// Using == as the operator is just for convenience.  Any binary operation
// that does not occur in the expressions could have been used.
//
ex to_nested_tuple(const lst& exprs, const lst& formulas)
{
    lst all;

    for (lst::const_iterator iter = exprs.begin(); iter != exprs.end(); ++iter) {
        all.append(*iter);
    }
    for (lst::const_iterator iter = formulas.begin(); iter != formulas.end(); ++iter) {
        all.append(*iter);
    }

    int n = all.nops();
    assert(n != 0);

    if (n == 1) {
        return all.op(0);
    }
    ex t = all.op(n-2) == all.op(n-1);
    for (int k = n - 3; k >= 0; --k) {
        t = all.op(k) == t;
    }
    return t;
}

void CloseContext(void)
{
 
	if (initialized) {

		for (lst::iterator i= MyHandleList.begin(); i != MyHandleList.end(); i++) {
			ThreadData* Tmp = *i;
			Tmp->thread_running = false;
				
		}
		Sleep(100);
		zmq_ctx_destroy (context);
		initialized = false;
	}
}

std::pair<ex,ex>  hyper_cube_den(lst pole_list,lst w_list, ex den)
{
try
  {
    using namespace lemon;
    GlpkLp lp;
    exhashmap<LpBase::Col> col_map;
    for(lst::const_iterator wi = w_list.begin();wi!=w_list.end();++wi)
      {
        col_map[*wi] = lp.addCol();
      }
    for(lst::const_iterator pit = pole_list.begin();pit!= pole_list.end();++pit)
      {
        ex tmp_expr = *pit;
        // cout<<"Expr: "<<*pit<<" subexprs:"<<endl;
        Lp::Expr constr_expr;
        for(lst::const_iterator wi = w_list.begin();wi!=w_list.end();++wi)
          {
            // cout<<*wi<<"  coeff "<<tmp_expr.coeff(*wi)<<endl;
            ex wi_coeff = tmp_expr.coeff(*wi);
            tmp_expr-=(*wi)*wi_coeff;
            if(is_a<numeric>(wi_coeff))
              {
                constr_expr+=ex_to<numeric>(wi_coeff).to_double()*col_map[*wi];
                //                cout<<ex_to<numeric>(wi_coeff).to_double()<<endl; 
              }
            else throw std::logic_error("Non numeric coefficient in pole term. ");
          }
        //constr_expr+=
        //cout<<"Ostatok "<<tmp_expr<<endl;
        if(is_a<numeric>(tmp_expr))
          lp.addRow(-ex_to<numeric>(tmp_expr).to_double(),constr_expr,Lp::INF);
        else throw std::logic_error(std::string("Lower bound is not a numeric"));
      }

    double l_bound,r_bound;
    cout<<"Hyper cube"<<endl;
    //   for(lst::const_iterator wi = w_list.begin();wi!=w_list.end();++wi)
      
        lp.min();
        lp.obj(col_map[*w_list.begin()]);
        // Solve the problem using the underlying LP solver
        lp.solve();
        if (lp.primalType() == Lp::OPTIMAL) 
          {
            cout<<*w_list.begin()<<" = ("<<lp.primal()<<",";
            l_bound = lp.primal();
          }
        else throw std::logic_error("Optimal solution not found.");
        lp.max();
        lp.obj(col_map[*w_list.begin()]);
        // Solve the problem using the underlying LP solver
        lp.solve();
        if (lp.primalType() == Lp::OPTIMAL) 
          {
            cout<<lp.primal()<<")"<<endl;
            r_bound = lp.primal();
          }
        else throw  std::logic_error("Optimal solution not found.");

        cout<<"Simplex : "<<endl;
        ex lbs,ubs;
        //        bool founds;
        //        tie(lbs,ubs,founds) = simplex_ara(pole_list,w_list,*w_list.begin());
        // cout<<"lbs: "<<lbs<<" ubs: "<<ubs<<endl;
        //  boost::mt19937 rng;       
        //        boost::uniform_real<> bounded_distribution(l_bound,r_bound);      // distribution that maps to 1..6
        // see random number distributions
        //  boost::variate_generator<boost::mt19937&, boost::uniform_real<> >
        //  die(rng, bounded_distribution);             // glues randomness with mapping
        //rng.seed(static_cast<unsigned char>(std::time(0)));

        
        //        double x   = die();                      // simulate rolling a die
        // new edition in center of interval no random
        //if((l_bound - trunc(l_bound)) < 10e-8) lbs = trunc(l_bound);
        //lbs = ginac_set_d(l_bound);
	lbs = l_bound;
        //if((r_bound - trunc(r_bound)) < 10e-8) ubs = trunc(r_bound);
	//ubs = ginac_set_d(r_bound);
	ubs = r_bound;
        //cout<<"ginac_set_d "<<ginac_set_d(l_bound)<<" "<<ginac_set_d(r_bound)<<endl;
        ex x_half = lbs*(1- pow(den,-1)) + ubs*pow(den,-1);
        //        double real_half = (r_bound + l_bound)/2.0;
        //        cout<<"Point "<<die()<<"   "<<x<<endl;
        return std::make_pair(*w_list.begin(),x_half);
  }catch(std::exception &p)
    {
      throw std::logic_error(std::string("In function \"hyper_cube_den\":\n |___> ")+p.what());
    }
}

/**
 *
 *  loop momentums,propagator expressions,
 *  invariants substitutions,propagator powers,number of loops
 *
 \param k_lst loop momentums list
 \param p_lst propagator expressions list
 \param subs_lst invariants substitutions list
 \param nu propagator powers list
 \param l number of loops 
   
 \return 
 *
 */
RoMB_loop_by_loop:: RoMB_loop_by_loop(
				      lst k_lst,
				      lst p_lst,
				      lst subs_lst,
                                      lst nu,
                                      bool subs_U
                                      )
{
  try
    {
      /* 
	 empty integral
      */
      MBintegral MBlbl_int(lst(),lst(),1);
      /* 
	 Full set of unused propagators, will change 
      */
      exlist input_prop_set;//( p_lst.begin(),p_lst.end());
      /* 
	 map for propagator powers
      */
      exmap prop_pow_map;
      for(lst::const_iterator Pit = p_lst.begin(); Pit != p_lst.end(); ++Pit)
	{
          input_prop_set.push_back(Pit->expand());
          prop_pow_map[Pit->expand()] = nu.op(std::distance(p_lst.begin(),Pit));
	}
	
      cout<<"INPSET: "<<input_prop_set<<endl;

      /* 
	 Iterate over momentums k1,k2,k3,etc.
      */
      unsigned int displacement_x = 0;
      unsigned int displacement_w = 0;
      for(lst::const_iterator kit = k_lst.begin(); kit != k_lst.end(); ++kit)
	{
          // Integral Normalization coefficient 
        //            MBlbl_int *= pow(I,k_lst.nops());
           MBlbl_int *= 1/tgamma(1+get_symbol("eps"));
          //MBlbl_int *= pow(Pi,2-get_symbol("eps"));
          //MBlbl_int *= exp(Euler*get_symbol("eps"));	  
          cout<<"PROP_POW_MAP "<<prop_pow_map<<endl;
	  /*
	    temporary set of propagators, with all momentum,except deleted
	  */
	  exlist tmp_p_lst(input_prop_set.begin(), input_prop_set.end()); 
	  /*
	    temporary set of propagators, with KIT momentum
	  */
	  lst P_with_k_lst;
	  BOOST_FOREACH(ex prop_tmp, tmp_p_lst)
	    {
	      if(prop_tmp.has(*kit))
		{
		  P_with_k_lst.append(prop_tmp);
		  input_prop_set.remove(prop_tmp);
		}
	    }
            
	  cout<< "Set wo k_i "<<input_prop_set<<endl;
	  cout<<" PWKlst "<<P_with_k_lst<<endl;
	  bool direct_formula_applied = false;
          // if only one term in PWKLST use well known formulas
          // [Smirnov A.1]
          if(!direct_formula_applied && (P_with_k_lst.nops() == 1))
            {
              ex pr_t = P_with_k_lst.op(0);
              ex nu_t = prop_pow_map[pr_t];
              exmap repls;
              BOOST_ASSERT_MSG(pr_t.match(-pow(*kit,2) + wild(2)),"ONE PROP");
              if(pr_t.match(-pow(*kit,2) + wild(2),repls))cout<<"repls: "<<repls<<endl;
              ex mass_tadpole = (tgamma(nu_t+get_symbol("eps")-2)/tgamma(nu_t)*pow(wild(2).subs(repls),-nu_t-get_symbol("eps")+2));
              cout<<mass_tadpole<<endl;
              MBlbl_int *= mass_tadpole;
              MBlbl_int.add_pole(nu_t+get_symbol("eps")-2);
              direct_formula_applied = true;
            }
          if(!direct_formula_applied && (P_with_k_lst.nops() == 2)) 
            {
              //TWO terms in PWK_LST, [Smirnov A.4]
              exmap repls_tad;
              if((P_with_k_lst.nops()==2) &&
        	 (( (P_with_k_lst.op(0).match(-pow(*kit,2))) && (P_with_k_lst.op(1).match(-pow(*kit,2)+wild())))||
                  ( (P_with_k_lst.op(1).match(-pow(*kit,2))) && (P_with_k_lst.op(0).match(-pow(*kit,2)+wild()))))
                 && !wild().has(*kit)
//.........这里部分代码省略.........