C++ ex::subs方法代码示例

// Solve
vector<double> Solve(GiNaC::ex Ex, GiNaC::symbol Sym) {
  vector<double> Roots;

  unsigned Degree = Ex.degree(Sym);
  auto Coeffs = GetCoeffs(Ex, Sym);
  // Bhaskara.
  if (Degree == 2) {
    GiNaC::ex A = Coeffs[2];
    GiNaC::ex B = Coeffs[1];
    GiNaC::ex C = Coeffs[0];

    GiNaC::ex Delta = B*B - 4 * A * C;
    // Guaranteed real roots.
    if (GiNaC::is_a<GiNaC::numeric>(Delta) &&
        !GiNaC::ex_to<GiNaC::numeric>(Delta).is_negative()) {
      GiNaC::ex Delta = GiNaC::sqrt(B*B - 4 * A * C).evalf();
      GiNaC::ex One = ((-B) + Delta)/(2*A);
      GiNaC::ex Two = ((-B) - Delta)/(2*A);
      if (GiNaC::is_a<GiNaC::numeric>(One))
        Roots.push_back(GiNaC::ex_to<GiNaC::numeric>(One.evalf()).to_double());
      if (GiNaC::is_a<GiNaC::numeric>(Two))
        Roots.push_back(GiNaC::ex_to<GiNaC::numeric>(Two.evalf()).to_double());
    }
  }
  // Cardano.
  else if (Degree == 3) {
    GiNaC::ex A = Coeffs[3];
    GiNaC::ex B = Coeffs[2];
    GiNaC::ex C = Coeffs[1];
    GiNaC::ex D = Coeffs[1];

    GiNaC::ex Delta0 = B*B - 3 * A * C;
    GiNaC::ex Delta1 = 2 * B*B*B - 9 * A * B * C + 27 * A * A * D;
    GiNaC::ex CD = Delta1 + GiNaC::sqrt(Delta1 * Delta1 - 4 * GiNaC::pow(Delta0, 3));
    CD = CD/2;
    CD = GiNaC::pow(CD, GiNaC::numeric(1)/3);

    GiNaC::symbol U("u");
    GiNaC::ex Var = GiNaC::numeric(-1)/(3 * A) * (B + U * CD + Delta0/(U * CD));
    GiNaC::ex One   = Var.subs(U == 1);
    GiNaC::ex Two   = Var.subs(U == ((-1 + GiNaC::sqrt(GiNaC::numeric(-3)))/2));
    GiNaC::ex Three = Var.subs(U == ((-1 - GiNaC::sqrt(GiNaC::numeric(-3)))/2));
    if (GiNaC::is_a<GiNaC::numeric>(One))
      Roots.push_back(GiNaC::ex_to<GiNaC::numeric>(One.evalf()).to_double());
    if (GiNaC::is_a<GiNaC::numeric>(Two))
      Roots.push_back(GiNaC::ex_to<GiNaC::numeric>(Two.evalf()).to_double());
    if (GiNaC::is_a<GiNaC::numeric>(Three))
      Roots.push_back(GiNaC::ex_to<GiNaC::numeric>(Three.evalf()).to_double());
  }

  return Roots;
}

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

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

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

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

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