C++ ginac::ex类代码示例

    // treat powers specially, because it's better to unroll them and
    // avoid a function call if we have an exact integral power
    std::string cpp_cse::print_power(const GiNaC::ex& expr, bool use_count)
      {
        std::ostringstream out;
        size_t n = expr.nops();

        if(n != 2)
          {
            error_context err_ctx = this->data_payload.make_error_context();
            err_ctx.error(ERROR_CSE_POWER_ARGUMENTS);
            out << "std::pow(" << this->print_operands(expr, ",", use_count) << ")";
            return std::string{};
          }

        // set up aliases for base and exponent expressions
        const auto& base_expr = expr.op(0);
        const auto& exponent_expr = expr.op(1);
    
        // perform use-counting on exponent, which is necessary since its values may have been
        // used elsewhere in the CSE tree, even if it is an integer that we will unroll and not use explicitly
        std::string exponent;
        if(use_count) exponent = this->get_symbol_with_use_count(exponent_expr);
        else exponent = this->get_symbol_without_use_count(exponent_expr);

        if(GiNaC::is_a<GiNaC::numeric>(exponent_expr))
          {
            const auto& exp_numeric = GiNaC::ex_to<GiNaC::numeric>(exponent_expr);
        
            std::string base;
            if(use_count) base = this->get_symbol_with_use_count(expr.op(0));
            else base = this->get_symbol_without_use_count(expr.op(0));

            if(GiNaC::is_integer(exp_numeric))
              {
                if(GiNaC::is_nonneg_integer(exp_numeric))
                  {
                    if(exp_numeric.to_int() == 0) out << "1.0";
                    else if(exp_numeric.to_int() == 1) out << unwrap_power(base, base_expr, 1);
                    else if(exp_numeric.to_int() == 2) out << "(" << unwrap_power(base, base_expr, 2) << ")";
                    else if(exp_numeric.to_int() == 3) out << "(" << unwrap_power(base, base_expr, 3) << ")";
                    else if(exp_numeric.to_int() == 4) out << "(" << unwrap_power(base, base_expr, 4) << ")";
                    else out << "std::pow(" << base << "," << exp_numeric.to_int() << ")";
                  }
                else  // negative integer
                  {
                    if(exp_numeric.to_int() == -0) out << "1.0";
                    else if(exp_numeric.to_int() == -1) out << "1.0/" << unwrap_power(base, base_expr, 1);
                    else if(exp_numeric.to_int() == -2) out << "1.0/" << "(" << unwrap_power(base, base_expr, 2) << ")";
                    else if(exp_numeric.to_int() == -3) out << "1.0/" << "(" << unwrap_power(base, base_expr, 3) << ")";
                    else if(exp_numeric.to_int() == -4) out << "1.0/" << "(" << unwrap_power(base, base_expr, 4) << ")";
                    else out << "std::pow(" << base << "," << exp_numeric.to_int() << ")";
                  }
              }
            else  // not an integer
              {
                out << "std::pow(" << base << "," << exponent << ")";
              }
          }

        return(out.str());
      }

// GetSymbols
set<GiNaC::symbol, GiNaC::less> GetSymbols(GiNaC::ex Ex) {
  set<GiNaC::symbol, GiNaC::less> Syms;
  for (auto EI = Ex.preorder_begin(), E = Ex.preorder_end(); EI != E; ++EI)
    if (GiNaC::is_a<GiNaC::symbol>(*EI))
      Syms.insert(GiNaC::ex_to<GiNaC::symbol>(*EI));

  return Syms;
}

// GetCoeffs
vector<GiNaC::ex> GetCoeffs(GiNaC::ex Ex, GiNaC::symbol Sym) {
  assert(Ex.is_polynomial(Sym));
  vector<GiNaC::ex> Coeffs;
  unsigned Degree = Ex.degree(Sym);
  for (unsigned Idx = 0; Idx <= Degree; ++Idx)
    Coeffs.push_back(Ex.coeff(Sym, Idx));
  return Coeffs;
}

Expr Expr::subs(std::vector<std::pair<Expr, Expr> > Subs) { 
  GiNaC::ex Ex = Expr_;
  for (auto& E : Subs) {
    //dbgs() << "Replacing " << E.first.getExpr() << " with " << E.second.getExpr() << " in " << Ex;
    Ex = Ex.subs(E.first.getExpr() == E.second.getExpr());
    //dbgs() << " giving " << Ex << "\n";
  }
  return Expr(Ex);
}

 const flattened_tensor& CovariantRiemannB3Cache::get()
   {
     if(this->B3) return *this->B3;
 
     auto args = res.generate_cache_arguments(printer);
 
     this->B3 = std::make_unique<flattened_tensor>(res.fl.get_flattened_size<field_index>(RESOURCE_INDICES::RIEMANN_B3_INDICES));
 
     const auto max = res.share.get_max_field_index(variance::covariant);
     const auto max_l = res.share.get_max_field_index(variance::contravariant);
 
     SubstitutionMapCache subs_cache(res, printer);
     DerivativeSymbolsCache deriv_cache(res, res.share, printer);
 
     for(field_index i = field_index(0, variance::covariant); i < max; ++i)
       {
         for(field_index j = field_index(0, variance::covariant); j < max; ++j)
           {
             for(field_index k = field_index(0, variance::covariant); k < max; ++k)
               {
                 unsigned int index = res.fl.flatten(i,j,k);
 
                 GiNaC::ex subs_expr = 0;
 
                 if(!res.cache.query(expression_item_types::Riemann_B3_item, index, args, subs_expr))
                   {
                     timing_instrument timer(res.compute_timer);
 
                     auto& deriv_syms = deriv_cache.get();
                     GiNaC::ex expr = 0;
 
                     for(field_index l = field_index(0, variance::contravariant); l < max_l; ++l)
                       {
                         auto Rie_ijk = (*res.Rie_T)(static_cast<unsigned int>(k), static_cast<unsigned int>(i),
                                                     static_cast<unsigned int>(j), static_cast<unsigned int>(l));
                         auto Rie_ikj = (*res.Rie_T)(static_cast<unsigned int>(k), static_cast<unsigned int>(j),
                                                     static_cast<unsigned int>(i), static_cast<unsigned int>(l));
                         
                         auto Rie_sym = (Rie_ijk + Rie_ikj) / 2;
                         
                         expr += Rie_sym * deriv_syms[res.fl.flatten(l)];
                       }
 
                     // get substitution map
                     GiNaC::exmap& subs_map = subs_cache.get();
                     subs_expr = expr.subs(subs_map, GiNaC::subs_options::no_pattern);
                     
                     res.cache.store(expression_item_types::Riemann_B3_item, index, args, subs_expr);
                   }
                 
                 (*this->B3)[index] = subs_expr;
               }
           }
       }
 
     return *this->B3;
   }

int Expr::getMaxDegree() const {
  int Max = 0;
  for (auto It = Expr_.preorder_begin(), E = Expr_.preorder_end();
       It != E; ++It) {
    GiNaC::ex Ex = *It;
    if (GiNaC::is_a<GiNaC::power>(Ex))
      Max = std::max(Max, GiNaC::ex_to<GiNaC::numeric>(Ex.op(1)).to_int());
    else if (GiNaC::is_a<GiNaC::symbol>(Ex))
      Max = std::max(Max, 1);
  }
  return Max;
}

  MatrixWrapper::SymmetricMatrix
  NonLinearAnalyticConditionalGaussian_Ginac::CovarianceGet() const
  {
    if (cond_size!=0)
      {
	MatrixWrapper::ColumnVector u_num   (u_size);
	MatrixWrapper::ColumnVector x_num   (x_size);
	GiNaC::ex substitute (func_size);
	MatrixWrapper::Matrix D             (func_size,cond_size);


	u_num = ConditionalArgumentGet(1);
	x_num = ConditionalArgumentGet(0);

	for (unsigned int i=0; i<cond_size; i++)
	  {
	    // temp variable to substitute in
	    substitute = dfunc_dcond[i];

	    // substitute all u_sym with u_num
	    for (unsigned int j=0; j<u_size; j++)
	      substitute = substitute.subs( u_sym[j]==u_num(j+1) );

	    // substitute all x_sym with x_num
	    for (unsigned int j=0; j<x_size; j++)
	      substitute = substitute.subs( x_sym[j]==x_num(j+1) );

	    // convert substitute back to matrix
	    GiNaC::matrix substitute_matrix = GiNaC::ex_to<GiNaC::matrix>(substitute);

	    // build matrix D
	    for (unsigned int j=0; j<func_size; j++)
	      D(j+1,i+1) = GiNaC::ex_to<GiNaC::numeric>( substitute_matrix(j,0).evalf() ).to_double();
	  }
	//cout << "D: " << D << endl;
	//cout << "CondCov:\n" << (Matrix)cond_covariance << endl;
	MatrixWrapper::Matrix temp = D * (MatrixWrapper::Matrix)AdditiveNoiseSigmaGet() * D.transpose();

	// convert func_covariance_matrix to symmetric matrix
	MatrixWrapper::SymmetricMatrix additiveNoise(temp.rows());
	temp.convertToSymmetricMatrix(additiveNoise);
	return additiveNoise;

      }
    else
      {
	return AdditiveNoiseSigmaGet();
      }
  }

GiNaC::ex operator_from_graph(KontsevichGraph graph, PoissonStructure poisson, std::vector<GiNaC::ex> arguments)
{
    // TODO: check if graph.external() == arguments.size()
    GiNaC::ex result = 0;
    map_operator_coefficients_from_graph(graph, poisson, [&result, &poisson, &arguments](std::vector< std::multiset<size_t> > derivatives, GiNaC::ex summand) {
        GiNaC::ex tail = 1;
        for (size_t m = 0; m != arguments.size(); ++m)
        {
            GiNaC::ex factor = arguments[m];
            for (size_t idx : derivatives[m])
                factor = factor.diff(poisson.coordinates[idx]);
            tail *= factor;
        }
        result += summand * tail;
    });
    return result;
}

  MatrixWrapper::ColumnVector
  NonLinearAnalyticConditionalGaussian_Ginac::ExpectedValueGet() const
  {
    MatrixWrapper::ColumnVector u_num   (u_size);
    MatrixWrapper::ColumnVector x_num   (x_size);
    MatrixWrapper::ColumnVector func_num(func_size);
    GiNaC::ex substitute (func_size);
    MatrixWrapper::ColumnVector expected(func_size);

    u_num = ConditionalArgumentGet(1);
    x_num = ConditionalArgumentGet(0);

    // use Mu of additive noise
    if (cond_size!=0)
      for (unsigned int i=0; i<u_size; i++)
	for (unsigned int j=0; j<cond_size; j++)
	  if (u_sym[i] == cond_sym[j])
	      u_num(i+1) += (this->AdditiveNoiseMuGet())(j+1);


    // evaluate func
    for (unsigned int i=0; i<func_size; i++)
      {
	// temp variable to substitute in
	substitute = func_sym(i,0);

	// substitute all u_sym with u_num
	for (unsigned int j=0; j<u_size; j++)
	  substitute = substitute.subs( u_sym[j]==u_num(j+1) );

	// substitute all x_sym with x_num
	for (unsigned int j=0; j<x_size; j++)
	  substitute = substitute.subs( x_sym[j]==x_num(j+1) );

	// build matrix func_num
	func_num(i+1) = GiNaC::ex_to<GiNaC::numeric>( substitute.evalf() ).to_double();
      }
    expected = func_num;

    if (cond_size==0)
      expected += AdditiveNoiseMuGet();

    return expected;
  }

GiNaC::lst VectorField::FindVarsInEx(const GiNaC::ex &e)
    {
    unsigned k;
    GiNaC::lst vlist;
    for (k = 0; k < varname_list.nops(); ++k)
        {
        if (e.has(varname_list[k]))
            vlist.append(varname_list[k]);
        }
    return vlist;
    }

// NegativeOrdinate
// Returns the ranges for which the ordinate is negative.
vector<pair<GiNaC::ex, GiNaC::ex> >
    NegativeOrdinate(vector<double> Vec, GiNaC::symbol Sym, GiNaC::ex Ex) {
  vector<pair<GiNaC::ex, GiNaC::ex> > Negs;

  // No real roots - it's either all positive or all negative.
  if (Vec.empty()) {
    GiNaC::ex Subs = Ex.subs(Sym == 0);
    if (GiNaC::is_a<GiNaC::numeric>(Subs) &&
        GiNaC::ex_to<GiNaC::numeric>(Subs).is_negative())
        Negs.push_back(make_pair(GiNaC::inf(-1), GiNaC::inf(1)));
    return Negs;
  }

  sort(Vec.begin(), Vec.end());

  GiNaC::ex FirstSubs = Ex.subs(Sym == GiNaC::numeric(Round(Vec[0] - 1.0f)));
  // Negative at -inf to the first root minus one.
  if (GiNaC::is_a<GiNaC::numeric>(FirstSubs) &&
      GiNaC::ex_to<GiNaC::numeric>(FirstSubs).is_negative())
    Negs.push_back(make_pair(GiNaC::inf(-1), GiNaC::numeric(Vec[0] - 1.0f)));

  for (unsigned Idx = 1; Idx < Vec.size(); ++Idx) {
    long int E = Round(Vec[Idx]);
    // Check if it's negative to the right.
    GiNaC::ex Subs = Ex.subs(Sym == GiNaC::numeric(E + 1));
    if (GiNaC::is_a<GiNaC::numeric>(Subs))
      if (GiNaC::ex_to<GiNaC::numeric>(Subs).is_negative()) {
        // Last root, it's negative all the way to +inf.
        if (Idx == Vec.size() - 1)
          Negs.push_back(make_pair(GiNaC::numeric(E + 1), GiNaC::inf(1)));
        // Not the last root, it's negative from this root to the next.
        else
          Negs.push_back(make_pair(GiNaC::numeric(E + 1),
                                   GiNaC::numeric(Round(Vec[Idx + 1] - 1))));
      }
  }

  return Negs;
}

bool TestExpr::runOnModule(Module&) {
  Expr A("a");
  Expr B("b");

  assert(A + B == B + A);
  assert(A * B == B * A);
  assert(A - B != B - A);
  assert(A / B != B / A);

  //map<GiNaC::symbol, int> GetMaxExps(GiNaC::ex Ex); 
  GiNaC::symbol GA("a");
  GiNaC::symbol GB("b");
  GiNaC::ex EZ = 5 * GA * GA * GA + 4 * GB + 1;

  assert(EZ.is_polynomial(GA));
  assert(EZ.degree(GA) == 3);
  assert(EZ.coeff(GA, 3) == GiNaC::numeric(5));
  assert(EZ.coeff(GA, 2) == GiNaC::numeric(0));
  assert(EZ.coeff(GA, 1) == GiNaC::numeric(0));
  assert(EZ.coeff(GA, 0) == 4 * GB + 1);

  assert(EZ.is_polynomial(GB));
  assert(EZ.degree(GB) == 1);
  assert(EZ.coeff(GB, 1) == GiNaC::numeric(4));
  assert(EZ.coeff(GB, 0) == 5 * GA * GA * GA + 1);

  auto SEZ = GetSymbols(EZ);
  assert(SEZ.size() == 2);
  assert(SEZ.count(GA) && SEZ.count(GB));

  auto CEZGA = GetCoeffs(EZ, GA);
  assert(CEZGA.size() == 4);
  assert(CEZGA[3] == GiNaC::numeric(5));
  assert(CEZGA[2] == GiNaC::numeric(0));
  assert(CEZGA[1] == GiNaC::numeric(0));
  assert(CEZGA[0] == 4 * GB + 1);

  auto CEZGB = GetCoeffs(EZ, GB);
  assert(CEZGB.size() == 2);
  assert(CEZGB[1] == GiNaC::numeric(4));
  assert(CEZGB[0] == 5 * GA * GA * GA + 1);

  auto EF = GA * GA - 1;
  auto EFSolve = Solve(EF, GA);
  assert(EFSolve.size() == 2);
  assert(Round(EFSolve[0]) == -1 || Round(EFSolve[1]) == -1);
  assert(Round(EFSolve[0]) == 1  || Round(EFSolve[1]) == 1);
  auto EFNegs = NegativeOrdinate(EFSolve, GA, EF); 
  assert(EFNegs.empty());

  EF = GA * GA - 2 * GA + 1;
  EFSolve = Solve(EF, GA);
  assert(EFSolve.size() == 2);
  assert(Round(EFSolve[0]) == 1 && Round(EFSolve[1]) == 1);
  EFNegs = NegativeOrdinate(EFSolve, GA, EF); 
  assert(EFNegs.empty());

  EF = (-1) * GA * GA + GA - 1;
  EFSolve = Solve(EF, GA);
  assert(EFSolve.size() == 0);
  EFNegs = NegativeOrdinate(EFSolve, GA, EF); 
  assert(EFNegs.size() == 1);
  assert(EFNegs[0].first.is_equal(GiNaC::inf(-1)));
  assert(EFNegs[0].second.is_equal(GiNaC::inf(1)));

  return false;
}

 bool operator()(GiNaC::ex const& l, GiNaC::ex const& r) const {
   return l.compare(r) < 0;
 }

double Attribute::EvalCompiledExpression (double const val, std::string const attrib ) {

//cout << GetPrototype()->GetName() << " ??  at pointer num " << m_cur_fp << " -> compiled = " << m_compiled.at(m_cur_fp) << endl;
 	if (!m_compiled.at(m_cur_fp)) {
 		//substitute all attributes with numbers in GiNaC expression, except the attribute
 		//which serves as the free parameter for runtime compilation
 		GiNaC::lst symlist;
 		GiNaC::lst numlist;
 		for (unsigned int i=0; i<m_subjects.size() ; i++) {
 			Attribute* a = m_subjects.at(i);
 			if (a->GetName() == attrib) continue;
 	        symlist.append( get_symbol(a->GetSymbol()) );
 			if (a->GetTypeID()==typeid(  double*).name()) { numlist.append(a->GetMember  <double>() ); continue; }
 			if (a->GetTypeID()==typeid(     int*).name()) { numlist.append(a->GetMember     <int>() ); continue; }
 			if (a->GetTypeID()==typeid(    long*).name()) { numlist.append(a->GetMember    <long>() ); continue; }
 			if (a->GetTypeID()==typeid(unsigned*).name()) { numlist.append(a->GetMember<unsigned>() ); continue; }
 			if (a->GetTypeID()==typeid(    bool*).name()) { numlist.append(a->GetMember    <bool>() ); continue; }
 		}

		GiNaC::ex e = GiNaC::evalf((symlist.nops()==0)?m_expression:m_expression.subs(symlist,numlist));

		//add function pointers
		m_fp.push_back(NULL);
		m_fpi.push_back(NULL);
		//compile the GiNaC expression
		try {
			//fairly easy for real valued expressions
			if (!m_complex) {
				compile_ex(e, get_symbol(GetPrototype()->GetAttribute(attrib)->GetSymbol()), m_fp.at(m_num_fp));
			}
			//more work to do, since GiNaC::realsymbol does not behave as expected (and it is therefore not used at all)
			else {
				stringstream se; se << e; std::string formula = se.str();
				std::string sym  = GetPrototype()->GetAttribute(attrib)->GetSymbol();
				std::string asym = "abs(VarForEvalCompiledExpression)";
				Prototype::ReplaceString(formula,sym,asym);
				GiNaC::lst symlist;
				symlist.append( get_symbol("VarForEvalCompiledExpression") );
				GiNaC::ex ea = GiNaC::ex(formula,symlist);
				symlist.remove_all();
				symlist.append( get_symbol(sym) );

				GiNaC::ex ear = ea.real_part();
				stringstream ser;  ser << ear; formula = ser.str();
				if ( Prototype::ReplaceString(formula,asym,sym) ) {
					ear  = GiNaC::ex(formula,symlist);
					compile_ex(ear, get_symbol(GetPrototype()->GetAttribute(attrib)->GetSymbol()), m_fp.at(m_num_fp));
				}

				GiNaC::ex eai = ea.imag_part();
				stringstream sei;  sei << eai; formula = sei.str();
				if ( Prototype::ReplaceString(formula,asym,sym) ) {
					eai  = GiNaC::ex(formula,symlist);
					compile_ex(eai, get_symbol(GetPrototype()->GetAttribute(attrib)->GetSymbol()), m_fpi.at(m_num_fp));
				}
			}
 			//cout << " compiling expression " << e << " of attribute " << GetName() << " in module " << GetPrototype()->GetName() << endl;
 		 	m_num_fp++;
		}
	 	catch (exception &p) {
 			cout << " Warning: attribute " << GetName() << " of module " << GetPrototype()->GetName() << endl << endl
 				 << " function Attribute::EvalCompiledExpression" << endl
 				 << " No external runtime compiler available: " << p.what() << endl
				 << " Falling back to (slow) analytic evaluation!" << endl << endl
				 << " Hint: if you have a shell and gcc on your system, create the one-liner " << endl << endl
				 << "    #!/bin/sh" << endl
				 << "    gcc -x c -fPIC -shared -o $1.so $1" << endl << endl
 		         << " name it \"ginac-excompiler\", and put it somewhere in your search path." << endl << endl;
	 		m_ginac_excomp = false;
	 	}
		m_compiled.at(m_cur_fp) = true; //even if compilation failed, as we don't have to try a second time!
 	}

	//if compilation failed, invoke slow analytic evaluation
	if (!m_ginac_excomp ) {
		*((double*) GetPrototype()->GetAttribute(attrib)-> GetAddress()) = val;
 		EvalExpression();
 		return *((double*) GetAddress());
 	}

	//invoke fast runtime compiled routines
 	if (m_fpi.at(m_cur_fp) != NULL ) m_imaginary = m_fpi.at(m_cur_fp)(val);
	if ( m_fp.at(m_cur_fp) != NULL ) return m_fp.at(m_cur_fp)(val);

	return 0.0;

}

double eval_at(const math::ex &expr, double x_val, double y_val)
{
    math::ex temp = expr.subs(x == x_val).subs(y == y_val).evalf();
    return math::ex_to<math::numeric>( temp ).to_double();
}