C++ std::isfinite方法代码示例

floating_t Pyramid::calculate_link(const Node& node, const Node& desc,
                                   spectral_distance_meassure_t spec_dst) const
{
    using std::exp;
    using std::isfinite;
    using std::isinf;
    if (link(node, desc) == 0.l
        /* || 1./util::condition(value2(desc)) < 1e-6*/) {
        return 0.l;
    }
    const auto beta_spec = -0.01;
    const auto beta_vari = -0.01;
    const auto beta_spat = -0.01;

    auto spec_term = exp(beta_spec * spec_dst(node, desc));
    auto spat_term = exp(beta_spat * spatial_distance(node, desc));
    auto vari_term = exp(beta_vari * var(desc));

    auto result = spec_term * vari_term * spat_term;

    assert(isfinite(spec_term));
    assert(isfinite(vari_term));
    assert(isfinite(spat_term));
    assert(isfinite(result));

    return result;
}

 EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a, const Scalar& b) const {
   using std::isfinite;
   if (isfinite(a) && isfinite(b)) return a + b;
   if (isfinite(a)) return a;
   if (isfinite(b)) return b;
   return a + b;
 }

/**
 * \brief TODO
 */
static double
trunc_estim_param_bisec(const double mean, const double tol){

  double lambda_low = mean-1;
  double lambda_high = mean;
  double lambda_mid = mean - 0.5;

  double diff = std::numeric_limits<double>::max();
  double prev_val = std::numeric_limits<double>::max();

  while(movement(lambda_high, lambda_low) > tol &&
         diff > tol){

    lambda_mid = (lambda_low + lambda_high)/2;

    const double mid_val = lambda_score_funct(mean, lambda_mid);

    if (mid_val < 0) lambda_low = lambda_mid;
    else lambda_high = lambda_mid;

    diff = fabs((prev_val - mid_val)/std::max(mid_val, prev_val));

    prev_val = mid_val;
  }

  if (!(isfinite(lambda_mid))) {
    stringstream ss;
    ss << "Zero-truncated Poisson parameter estimation failed. Reason: "
       << "got non-finite value for lambda";
    throw DistributionException(ss.str());
  }

  return lambda_mid;
}

    T invoke(I& integrand, P const& point)
    {
        using std::isfinite;

        hep::projector<T> projector(this, point);

        // call the integrand function with the supplied point. Distributions
        // are generated here
        T value = integrand.function()(point, projector);

        if (value != T())
        {
            value *= point.weight();

            if (isfinite(value))
            {
                accumulate(sums_[0], sums_[1], compensations_[0], value);
                ++finite_calls_[0];
            }
            else
            {
                value = T();
            }

            ++non_zero_calls_[0];
        }

        return value;
    }

int luaV_tostring (lua_State *L, StkId obj) {
  if (!ttisnumber(obj))
    return 0;
  else {
    char s[LUAI_MAXNUMBER2STR];
    lua_Number n = nvalue(obj);
    // SPRING -- synced safety change
    //        -- need a custom number formatter?
    if (isfinite(n)) {
      lua_number2str(s, n);
    }
    else {
      if (isnan(n)) {
        strcpy(s, "nan");
      }
      else {
        const int inf_type = isinf(n);
        if (inf_type == 1) {
          strcpy(s, "+inf");
        } else if (inf_type == -1) {
          strcpy(s, "-inf");
        } else {
          strcpy(s, "weird_number");
        }
      }
    } 
    setsvalue2s(L, obj, luaS_new(L, s));
    return 1;
  }
}

    void add_to_2d_distribution(std::size_t index, T x, T y, T value)
    {
        using std::isfinite;

        if (!isfinite(value))
        {
            return;
        }

        // TODO: index might be larger than the than allowed; throw?
        auto const parameters = parameters_.at(index);

        T const shifted_x = x - parameters.x_min();

        if (shifted_x < T())
        {
            // point is outside the binning range
            return;
        }

        T const shifted_y = y - parameters.y_min();

        if (shifted_y < T())
        {
            // point is outside the binning range
            return;
        }

        std::size_t const bin_x = shifted_x / parameters.bin_size_x();

        if (bin_x >= parameters.bins_x())
        {
            // point is right of the range that we are binning
            return;
        }

        std::size_t const bin_y = shifted_y / parameters.bin_size_y();

        if (bin_y >= parameters.bins_y())
        {
            return;
        }

        std::size_t const new_index = indices_.at(index) + 2 * (bin_y *
            parameters.bins_x() + bin_x);

        accumulate(
            sums_.at(new_index),
            sums_.at(new_index + 1),
            compensations_.at(new_index / 2),
            value
        );

        // FIXME: if this function is called more than once, the values are
        // incorrect
        ++non_zero_calls_.at(new_index / 2);
        ++finite_calls_.at(new_index / 2);
    }

void print_nonfinite(const FloatType& x, const Node& node, const Node& dst)
{
    using namespace std;
    if (!isfinite(x)) {
        cerr << endl
             << util::RED << "Non-finite weight " << x << "\tfrom " << node
             << "\tto " << dst << util::RESET << flush;
    }
}

bool TestExpDalpha(const MatrixXd& expDalpha) {
    for (int j = 0; j < expDalpha.cols(); j++) {
        for (int i = 0; i < expDalpha.rows(); i++) {
            if (isfinite(expDalpha(i, j)) == 0) {
                return false;
            }
        }
    }
    return true;
}

// check 3 term recurrence to avoid non-positive elements
// truncate if non-positive element found
static void
check_three_term_relation(vector<double> &a,
			  vector<double> &b){

  // first entry is zero! Abort
  if(a[0] <= 0.0){
    a.clear();
    b.clear();
  }

  for(size_t i = 0; i < b.size(); i++){
    if(b[i] <= 0.0 || !isfinite(b[i])
       || a[i + 1] <= 0.0 || !isfinite(a[i + 1])){
      b.resize(i);
      a.resize(i + 1);
      break;
    }
  }
}

void DPmatrix::compute_Pr_sum_all_paths()
{
  const int I = size1()-1;
  const int J = size2()-1;

  double total = 0.0;
  for(int state1=0;state1<nstates();state1++)
    total += (*this)(I,J,state1)*GQ(state1,endstate());

  Pr_total = pow(efloat_t(2.0),scale(I,J)) * total;
  assert(not isnan(log(Pr_total)) and isfinite(log(Pr_total)));
}

void ConvolvedSpikeTrain::UpdateExponentialVectors() {
  if (size > 0) {
    for (unsigned int i=0; i<size; ++i) {
      exp_pos[i] = exp((long double)(spikes[i]/tau));
      exp_neg[i] = exp((long double)(-spikes[i]/tau));
    }

    /* Check if any over/underflow occurred while calculating the
     * exponentials */
    if (tau!=0 && (exp_neg.front()==0 || !isfinite(exp_pos.back()))) {
      throw overflow_error("tau is too small compared to the spike times. Please use a larger value for tau, or shorter spike trains.");
    }
  }
}

bool
Modification::isFinite() const
{
    return  isfinite(m_transform.sx) &&
            isfinite(m_transform.shy) &&
            isfinite(m_transform.shx) &&
            isfinite(m_transform.sy) &&
            isfinite(m_transform.tx) &&
            isfinite(m_transform.ty);
}

floating_t Pyramid::calculate_var(const Node& node,
                                  spectral_distance_meassure_t spec_dst) const
{
    using accumulators::extract::sum;
    using common::descs;
    using accumulators::sum_acc;
    using std::isfinite;

    auto acc = sum_acc<floating_t>{};
    const auto node_norm = calculate_norm(node);

    for (const auto& desc : descs(node)) {
        if (link(node, desc) == 0.  // skip deterministic descs
            /* ||1./util::condition(value2(desc)) < 1e-6*/) {
            continue;
        }
        else {
            acc(spec_dst(node, desc) * calculate_weight(node, desc, node_norm));
        }
    }
    auto result = sum(acc);
    assert(isfinite(result));
    return result;
}

void Ip2_ViscElMat_ViscElMat_ViscElPhys::Calculate_ViscElMat_ViscElMat_ViscElPhys(const shared_ptr<Material>& b1, const shared_ptr<Material>& b2, const shared_ptr<Interaction>& interaction, shared_ptr<ViscElPhys> phys) {
	ViscElMat* mat1 = static_cast<ViscElMat*>(b1.get());
	ViscElMat* mat2 = static_cast<ViscElMat*>(b2.get());
	Real mass1 = 1.0;
	Real mass2 = 1.0;
	
	if ((isfinite(mat1->kn) and  not (isfinite(mat2->kn))) or
			(isfinite(mat2->kn) and  not (isfinite(mat1->kn))) or
			(isfinite(mat1->ks) and  not (isfinite(mat2->ks))) or
			(isfinite(mat2->ks) and  not (isfinite(mat1->ks))) or
			(isfinite(mat1->cn) and  not (isfinite(mat2->cn))) or
			(isfinite(mat2->cn) and  not (isfinite(mat1->cn))) or
			(isfinite(mat1->cs) and  not (isfinite(mat2->cs))) or
			(isfinite(mat2->cs) and  not (isfinite(mat1->cs))) or
			(isfinite(mat1->tc) and  not (isfinite(mat2->tc))) or
			(isfinite(mat2->tc) and  not (isfinite(mat1->tc))) or
			(isfinite(mat1->en) and  not (isfinite(mat2->en))) or
			(isfinite(mat2->en) and  not (isfinite(mat1->en))) or
			(isfinite(mat1->et) and  not (isfinite(mat2->et))) or
			(isfinite(mat2->et) and  not (isfinite(mat1->et)))) {
				throw runtime_error("Both materials should have the same defined set of variables e.g. tc, ks etc.!"); 
			}
			
	mass1 = Body::byId(interaction->getId1())->state->mass;
	mass2 = Body::byId(interaction->getId2())->state->mass;
	if (mass1 == 0.0 and mass2 > 0.0) {
		mass1 = mass2;
	} else if (mass2 == 0.0 and mass1 > 0.0) {
		mass2 = mass1;
	}
	
	// See [Pournin2001, just below equation (19)]
	const Real massR = mass1*mass2/(mass1+mass2);
	
	GenericSpheresContact* sphCont=YADE_CAST<GenericSpheresContact*>(interaction->geom.get());
	Real R1=sphCont->refR1>0?sphCont->refR1:sphCont->refR2;
	Real R2=sphCont->refR2>0?sphCont->refR2:sphCont->refR1;
	
	Real kn1 = 0.0; Real kn2 = 0.0;
	Real cn1 = 0.0; Real cn2 = 0.0;
	Real ks1 = 0.0; Real ks2 = 0.0;
	Real cs1 = 0.0; Real cs2 = 0.0;
	
	if (((isfinite(mat1->tc)) and (isfinite(mat1->en)) and (isfinite(mat1->et)))  or ((tc) and (en) and (et))) {
		//Set parameters according to [Pournin2001]
		
		const Real Tc = (tc) ? (*tc)(mat1->id,mat2->id) : (mat1->tc+mat2->tc)/2.0;
		const Real En = (en) ? (*en)(mat1->id,mat2->id) : (mat1->en+mat2->en)/2.0;
		const Real Et = (et) ? (*et)(mat1->id,mat2->id) : (mat1->et+mat2->et)/2.0;
    
    // Factor 2 at the end of each expression is necessary, because we calculate
    // individual kn1, kn2, ks1, ks2 etc., because kn1 = 2*kn, ks1 = 2*ks
    // http://www.mail-archive.com/[email protected]/msg08778.html
    kn1 = kn2 = 1/Tc/Tc * ( Mathr::PI*Mathr::PI + pow(log(En),2) )*massR*2;
    cn1 = cn2 = -2.0 /Tc * log(En)*massR*2;
    ks1 = ks2 = 2.0/7.0 /Tc/Tc * ( Mathr::PI*Mathr::PI + pow(log(Et),2) )*massR*2;
    cs1 = cs2 = -4.0/7.0 /Tc * log(Et)*massR*2;
    //           ^^^
    // It seems to be an error in [Pournin2001] (22) Eq.4, missing factor 2
    // Thanks to Dominik Boemer for pointing this out
    // http://www.mail-archive.com/[email protected]/msg08741.html

		if (std::abs(cn1) <= Mathr::ZERO_TOLERANCE ) cn1=0;
		if (std::abs(cn2) <= Mathr::ZERO_TOLERANCE ) cn2=0;
		if (std::abs(cs1) <= Mathr::ZERO_TOLERANCE ) cs1=0;
		if (std::abs(cs2) <= Mathr::ZERO_TOLERANCE ) cs2=0;
	} else if ((isfinite(mat1->kn)) and (isfinite(mat1->ks)) and (isfinite(mat1->cn)) and (isfinite(mat1->cs))) {
		//Set parameters explicitly
		kn1 = mat1->kn;
		kn2 = mat2->kn;
		ks1 = mat1->ks;
		ks2 = mat2->ks;
		cn1 = mat1->cn;
		cn2 = mat2->cn;
		cs1 = mat1->cs;
		cs2 = mat2->cs;
	} else {
		//Set parameters on the base of young modulus
		kn1 = 2*mat1->young*R1;
		kn2 = 2*mat2->young*R2;
		ks1 = kn1*mat1->poisson;
		ks2 = kn2*mat2->poisson;
		if ((isfinite(mat1->cn)) and (isfinite(mat1->cs))) {
			cn1 = mat1->cn;
			cn2 = mat2->cn;
			cs1 = mat1->cs;
			cs2 = mat2->cs;
		}
		else if( isfinite(mat1->en) and isfinite(mat1->et)) {
			const Real En = (en) ? (*en)(mat1->id,mat2->id) : (mat1->en+mat2->en)/2.0;
			cn1 = cn2 = 2.0*find_cn_from_en(En, massR,contactParameterCalculation(kn1,kn2),interaction);
			cs1 = cs2 = 0;
		}
	}
	
	const Real mR1 = mat1->mR;      const Real mR2 = mat2->mR; 
	const int mRtype1 = mat1->mRtype; const int mRtype2 = mat2->mRtype;
	
	
	phys->kn = contactParameterCalculation(kn1,kn2);
//.........这里部分代码省略.........