C++ lst类代码示例

void compile_ex(const lst& exprs, const lst& syms, FUNCP_CUBA& fp, const std::string filename)
{
	lst replacements;
	for (std::size_t count=0; count<syms.nops(); ++count) {
		std::ostringstream s;
		s << "a[" << count << "]";
		replacements.append(syms.op(count) == symbol(s.str()));
	}

	std::vector<ex> expr_with_cname;
	for (std::size_t count=0; count<exprs.nops(); ++count) {
		expr_with_cname.push_back(exprs.op(count).subs(replacements));
	}

	std::ofstream ofs;
	std::string unique_filename = filename;
	global_excompiler.create_src_file(unique_filename, ofs);

	ofs << "void compiled_ex(const int* an, const double a[], const int* fn, double f[])" << std::endl;
	ofs << "{" << std::endl;
	for (std::size_t count=0; count<exprs.nops(); ++count) {
		ofs << "f[" << count << "] = ";
		expr_with_cname[count].print(GiNaC::print_csrc_double(ofs));
		ofs << ";" << std::endl;
	}
	ofs << "}" << std::endl;

	ofs.close();

	global_excompiler.compile_src_file(unique_filename, filename.empty());
	// This is not standard compliant! ... no conversion between
	// pointer-to-functions and pointer-to-objects ...
	fp = (FUNCP_CUBA) global_excompiler.link_so_file(unique_filename+".so", filename.empty());
}

matrix getJacobian(const lst& vf, const lst& expr, const lst& params){
	int i, j, nv, np;
	
	nv = vf.nops();
	np = params.nops();
	
	matrix Jac = matrix(nv, np);
	
	for (i = 0; i < nv; ++i)
	{
		ex f;
		if (expr == 0)
			f = vf[i];
		else
			f= iterated_subs(vf[i],expr);
		
		for (j = 0; j < np; ++j){
			symbol v = ex_to<symbol>(params[j]);
            ex df = f.diff(v);
            if (df != 0)
				Jac(i,j) = df;
		}
	}
	return Jac;
}

//
// 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;
}

static void generate_lag_assignment(ofstream& fout, string name, const lst& lag,
        symbol& indvar, ex& history_formula)
{
    fout << "    if (" << indvar << " < " << lag.op(3) << ") {\n";
    fout << "        " << lag.op(0) << " = " << history_formula.subs(indvar == indvar - lag.op(3)) << endl;
    fout << "    }\n";
    fout << "    else {\n";
    fout << "        " << lag.op(0) << " = lagvalue(" << indvar - lag.op(3) <<  ", " <<  lag.op(1) << ")\n";
    fout << "    }\n";
}

/** Find all occurrences of a pattern. The found matches are appended to
 *  the "found" list. If the expression itself matches the pattern, the
 *  children are not further examined. This function returns true when any
 *  matches were found. */
bool ex::find(const ex & pattern, lst & found) const
{
	if (match(pattern)) {
		found.append(*this);
		found.sort();
		found.unique();
		return true;
	}
	bool any_found = false;
	for (size_t i=0; i<nops(); i++)
		if (op(i).find(pattern, found))
			any_found = true;
	return any_found;
}

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;
	}
}

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());
    }
}

ex iterated_subs(ex f, lst e)
{
    ex g = f;
    int n = e.nops();
    for (int i = n-1; i >= 0; --i) {
        g = g.subs(e[i]);
    }
    return g;
}

/** 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);
}

string p_to_deriv_string(lst vars, int p[])
{
    int n = vars.nops();
    ostringstream dname;
    for (int j = 0; j < n; ++j) {
        for (int k = 0; k < p[j]; ++k) {
            dname << "_" << vars[j];
        }
    }
    return dname.str();
}

//
// 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;
}

//
//  On return, lags is a GiNac::lst, where each element is a GiNac::lst
//  of length 4 containing {lagsym, variable_index + 1, var, lag_time}
//
void VectorField::convert_delay_to_lagvalue(ex& f, lst &lags)
{
    symbol t(IndependentVariable);
    exset dlist;
    f.find(delay(wild(1),wild(2)),dlist);
    for (exset::const_iterator iter = dlist.begin(); iter != dlist.end(); ++iter) {
        ex delayfunc = *iter;
        ex delayexpr = delayfunc.op(0);
        lst vars = FindVarsInEx(delayexpr);
        ex del = delayfunc.op(1);
        for (lst::const_iterator iter = vars.begin(); iter != vars.end(); ++iter) {
            ostringstream os;
            lst tmp;
            os << lags.nops() + 1;
            symbol lagsym("lag" + os.str());
            int vindex = FindVar(ex_to<symbol>(*iter));
            delayexpr = delayexpr.subs(*iter == lagsym);
            tmp = {lagsym, vindex + 1, *iter, del};
            lags.append(tmp);
        }
        f = f.subs(delayfunc == delayexpr);
    }
}

 GinacConstFunction::GinacConstFunction( const exvector& exs, const lst& params) 
     : Function( exs.size(), 0, params.nops()),
     exs_(exs), params_(params) {
         if(print__all){
         //std::cerr << "GinacConstFunction = ";					/**/
         //for(unsigned int i = 0; i < exs_.size(); i++ ){     /**/
         //    std::cerr << "[f("<<i<<") ="<<exs_[i] << "] " ; /**/
         //}													/**/
         //std::cerr << std::endl;								/**/
     	//}
     }
     
 void GinacConstFunction::eval(bool do_f, bool do_g, bool do_H){
     
     lst argsAndParams;
     lst::const_iterator params_it = params_.begin();
     
     for( unsigned int i = 0; params_it != params_.end(); params_it++, i++ ){
         argsAndParams.append( ex_to<symbol>(*params_it) == p(i));
     }
     //std::cerr << "NEXT FUNC"<< std::endl;           /**/
     
     if (do_f){
         for(unsigned int i = 0; i < exs_.size(); i++ ){
             f(i) = ex_to<numeric>(exs_[i].subs( argsAndParams )).to_double();
             //std::cerr << "f("<<i<<") ="<<f(i)<< std::endl; /**/
         }
     }
     
     if (do_g){
         clear_g();
     }
     
     if (do_H){
         clear_h();
     }
 }

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);
}