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

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

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

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

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

void generate_deriv(string lang, ofstream &fout, ofstream &pout, string name, int r,
                         lst vf, lst expreqn_list,
                         lst vars, lst params)
{
    if (r < 1) {
        cerr << "Error: generate_deriv called with order=" << r << ". order must be greater than 0.\n";
        exit(-1);
    }
    //
    // Print the function header
    //
    fout << endl;
    fout << "/*\n";
    if (r == 1) {
        fout << " *     deriv = Df(x)[v1]\n";
    }
    else if (r == 2) {
        fout << " *     deriv = D^2f(x)[v1,v2]\n";
    }
    else if (r == 3) {
        fout << " *     deriv = D^3f(x)[v1,v2,v3]\n";
    }
    else {
        fout << " *     deriv = D^" << r << "f(x)[v1,...,v" << r << "]\n";
    }
    fout << " */\n";
    if (lang == "c") {
        fout << "void " << name << "_diff" << r << "(double deriv[], double x_[], double p_[],";
        pout << "void " << name << "_diff" << r << "(double deriv[], double x_[], double p_[],";
    }
    else {
        fout << "function " << name << "_diff" << r << "(x_, p_,";
    }
    for (int k = 0; k < r; ++k) {
        if (lang == "c") {
            fout << "double v" << k+1 << "_[]";
            pout << "double v" << k+1 << "_[]";
        }
        else {
            fout << "v" << k+1 << "_";
        }
        if (k < r-1) {
            fout << ",";
            if (lang == "c") {
                pout << ",";
            }
        }
    }
    fout << ")\n";
    if (lang == "c") {
        pout << ");\n";
    }
    fout << "{\n";
    if (lang == "c") {
        CDeclare(fout,"double",vars);
        CDeclare(fout,"double",params);
        fout << endl;
    }
    const char* pre = "    ";
    if (lang == "javascript") {
        pre = "    var ";
    }
    GetFromVector(fout, pre, vars, "=", "x_", "[]", 0, ";");
    fout << endl;
    GetFromVector(fout,pre,params, "=", "p_", "[]", 0, ";");
    fout << endl;
/*
    int na = exprnames.nops();
    for (int i = 0; i < na; ++i) {
        fout << "    " << exprnames[i] << " = " << exprformulas[i] << ";" << endl;
    }
    if (na > 0) {
        fout << endl;
    }
*/

    if (lang == "javascript") {
        fout << "    var deriv = [];\n";
    }
    int n = vars.nops();
    int *p = new int[n];
    for (int i = 0; i < n; ++i) {
        if (lang == "c") {
            fout << "{\n";
        }
        // ex f = vf[i];
        ex f = iterated_subs(vf[i],expreqn_list);
        map<string,ex> derivatives;
        // Initialize p to [r,0,0,...,0]
        p[0] = r;
        for (int k = 1; k < n; ++k) {
            p[k] = 0;
        }
        fout << "    /*\n";
        fout << "     *  Partial derivatives of vf[" << i << "]. \n";
        fout << "     *  Any derivative not listed here is zero.\n";
        fout << "     */\n";
        do {
            // Compute the partial derivative
            ex df = f;
//.........这里部分代码省略.........