C++ VectorXd::sum方法代码示例

void parameters::Likelihood(const VectorXd & eff){
	 VectorXd loc;
	 loc.resize(m_proba.rows());
	 loc=(m_proba.rowwise().sum().array().log());
	 m_loglikelihood=eff.transpose()*loc;
	 m_bic=m_loglikelihood - 0.5*m_nbparam*log(eff.sum());
}

double unimcd_in(
		const VectorXd& m_resd,
		const int& h
	){
	const int n1=m_resd.size(),len=n1-h+1;
	double initmean=0.0,initcov=0.0,sumw=0.0;
	int minone;
	if(h==n1){
		initmean=m_resd.sum()/(double)h;
		initcov=(m_resd.array()-initmean).abs2().sum()/(double)(h-1);
		return(sqrt(initcov));
	}
	VectorXd y=m_resd;
	VectorXd ay(len);
	VectorXd ay2(len);
	VectorXd sq(len);
	VectorXd y2(n1);

	std::sort(y.data(),y.data()+y.size());
	ay(0)=y.head(h).sum();	
	for(int samp=1;samp<len;samp++) ay(samp)=ay(samp-1)-y(samp-1)+y(samp+h-1);
	ay2=ay.array().square()/(double)h;
	y2=y.array().square();
	sq(0)=y2.head(h).sum()-ay2(0);
	for(int samp=1;samp<len;samp++) sq(samp)=sq(samp-1)-y2(samp-1)+y2(samp+h-1)-ay2(samp)+ay2(samp-1);
	initcov=sq.minCoeff(&minone);
	initcov/=(double)(h-1);
	initmean=ay(minone)/(double)h;
	return(initmean);
}

SpMat make_C(SpMat chol_K_inv,VectorXd h2s, SpMat ZtZ){
  SpMat Ki = chol_K_inv.transpose() * chol_K_inv;
  SpMat C = ZtZ;
  C /= (1.0 - h2s.sum());
  C += Ki;
  return C;
}

void LoadMatrices(MatrixXd &A, VectorXd &W, VectorXd &p, string weight_type) {

    A = MatrixXd::Zero(n_n, n_n);
    MatrixXd Apn = MatrixXd::Zero(n_p, n_n);
    MatrixXd abs_Apn = MatrixXd::Zero(n_p, n_n);
    VectorXd k = VectorXd::Zero(n_n);
    p = VectorXd::Zero(n_n);
    W = VectorXd::Zero(n_p);

    // Topological incidence matrix
    for (unsigned int j = 0; j < n_p; j++) {
        int idx_n1 = wds->agelemek.at(pipe_idx.at(j))->Get_Cspe_Index();
        int idx_n2 = wds->agelemek.at(pipe_idx.at(j))->Get_Cspv_Index();
        Apn(j, idx_n1)++;
        Apn(j, idx_n2)--;
        abs_Apn(j, idx_n1)++;
        abs_Apn(j, idx_n2)++;
    }

    A = Apn.transpose() * Apn;
    for (unsigned i = 0; i < n_n; i++) {
        A(i, i) = 0;
        for (unsigned int j = 0; j < n_n; j++) {
            if (A(i, j) != 0) {
                A(i, j) = 1.;
                // k(i)++;
            }
        }
    }

    if (0 == strcmp(weight_type.c_str(), "topology")) {
        for (unsigned int i = 0; i < W.size(); i++)
            W(i) = 1;
    }

    if (0 == strcmp(weight_type.c_str(), "dp")) {
        double dp, dp_max = -1., weight_min = 0.0001;
        for (unsigned int i = 0; i < W.size(); i++) {
            // dp = wds->agelemek.at(pipe_idx.at(i))->Get_dprop("mass_flow_rate");
            dp = wds->agelemek.at(pipe_idx.at(i))->Get_dprop("headloss");
            // dp = wds->agelemek.at(pipe_idx.at(i))->Get_dprop("length");
            W(i) = abs(dp);
            if (W(i) > dp_max)
                dp_max = W(i);
            //cout<<endl<<wds->agelemek.at(pipe_idx.at(i))->Get_nev()<<" dp="<<dp;
        }
        for (unsigned int i = 0; i < W.size(); i++) {
            W(i) /= dp_max;
            if (W(i) < weight_min)
                W(i) = weight_min;
            // cout<<endl<<wds->agelemek.at(pipe_idx.at(i))->Get_nev()<<" weight="<<W(i);
        }
    }

    // Final computations
    p = abs_Apn.transpose() * W;
    sumW = W.sum();

}

double entropy(const VectorXd &values) {
  MatrixXd fy = values.array() / values.sum();
  double invlog = 1 / std::log(2);
  double epsilon = std::numeric_limits<double>::epsilon();
  // h = sum(fy * log2(fy+epsilon))
  double h =  (fy.array() * (fy.array() + epsilon).log().array() * invlog).sum();
  return (-1) * h;
}

VectorXd Zscores2Post(VectorXd& Zs){
    VectorXd post(Zs.size());
    VectorXd Zsq  = Zs.array().square();
    for(int i = 0; i < Zsq.size(); i ++){
        VectorXd Ztemp = (Zsq.array() - Zsq[i])/2;
        VectorXd Zexp = Ztemp.array().exp();
        post[i] = 1/Zexp.sum();
    }
    return(post);
}

SpMat make_Sigma(std::vector<SpMat> ZKZts, VectorXd h2s,double tol){
  int n = ZKZts[0].rows();
  int h = h2s.size();
  MatrixXd R(n,n);
  R.setZero();
  for(int i = 0; i < h; i++){
    R += h2s[i] * ZKZts[i];
  }
  R.diagonal().array() += (1.0-h2s.sum());
  return R.sparseView(0,tol);
}

void PF::propagate(VectorXd &particles, VectorXd &weights, int t)
{
  for(int i = 0; i < N; i++) {
	double p = particles[i];
  	boost::normal_distribution<> f_rng(p, f_std); // sample from f = q
  	particles[i] = f_rng(rng); 
  	boost::math::normal g(p, g_std); 
        weights[i] *= pdf(g, y[t]);
  }
  weights /= weights.sum();  
}

// exponentiates and normalizes a vector
void expAndNormalize(VectorXd& v)
{
    if (v.size() == 0) return;

    double maxValue = v[0];
    for (int i = 1; i < v.size(); i++) {
        if (v[i] > maxValue)
	    maxValue = v[i];
    }

    v = (v.cwise() - maxValue).cwise().exp();
    double Z = v.sum();
    v /= Z;
}

int argrand(const VectorXd &v)
{
    double cutoff = v.sum() * rand() / (double)RAND_MAX;
    double cumSum = 0.0;
    for (int i = 0; i < v.size(); i++) {
        cumSum += v[i];
        if (cumSum >= cutoff) {
            return i;
        }
    }

    SVL_LOG(SVL_LOG_FATAL, "bug");
    return -1;
}

// If there are two draws and the first has two tiles and the second -- three tiles,
// then sizes = c(2,3) and array = c(m_{1t_1},m_{1t_2},m_{2t_1},m_{2t_2},m_{2t_3})
bool dlmcell(const string &filename, const VectorXd &sizes, const vector<double> &array) {
  bool error = false;
  if (array.size()!=sizes.sum()) { error = true; return(error); }
  ofstream out(filename.c_str(),ofstream::out);
  if (!out.is_open()) { error = true; return(error); }
  vector<double>::const_iterator it = array.begin();
  for ( int i = 0 ; i < sizes.size() ; i++ ) {
    for ( int j = 0 ; j < sizes(i) ; j++ ) {
      out << fixed << setprecision(6) << *it << " "; it++;
    }
    out << endl;
  }
  out.close( );
  return(error);
}

/* The Variational Bayes Expectation step for each group.
 *
 *  mutable: Group assignment probabilities, qZj
 *  returns: The complete-data (X,Z) free energy E[log p(X,Z)/q(Z)] for group j.
 *  throws: invalid_argument rethrown from other functions.
 */
template <class W, class C> double vbexpectation (
    const MatrixXd& Xj,         // Observations in group J
    const W& weights,           // Group Weight parameter distribution
    const vector<C>& clusters,  // Cluster parameter distributions
    MatrixXd& qZj,              // Observations to group mixture assignments
    const bool sparse           // Do sparse updates to groups
    )
{
  const int K  = clusters.size(),
            Nj = Xj.rows();

  // Get log marginal weight likelihoods
  const ArrayXd E_logZ = weights.Elogweight();

  // Initialise and set K = 1 defaults for cluster counts
  ArrayXi Kful = ArrayXi::Zero(1), Kemp = ArrayXi::Zero(0);

  // Find empty clusters if sparse
  if ( (sparse == false) && (K > 1) )
    Kful = ArrayXi::LinSpaced(Sequential, K, 0, K-1);
  else if (sparse == true)
    arrfind((weights.getNk() >= ZEROCUTOFF), Kful, Kemp);

  const int nKful = Kful.size(),
            nKemp = Kemp.size();

  // Find Expectations of log joint observation probs -- allow sparse evaluation
  MatrixXd logqZj(Nj, nKful);

  for (int k = 0; k < nKful; ++k)
    logqZj.col(k) = E_logZ(Kful(k)) + clusters[Kful(k)].Eloglike(Xj).array();

  // Log normalisation constant of log observation likelihoods
  const VectorXd logZzj = logsumexp(logqZj);

  // Make sure qZ is the right size, this is a nop if it is
  qZj.resize(Nj, K);

  // Normalise and Compute Responsibilities -- again allow sparse evaluation
  for (int k = 0; k < nKful; ++k)
    qZj.col(Kful(k)) = ((logqZj.col(k) - logZzj).array().exp()).matrix();

  // Empty Cluster Responsabilities
  for (int k = 0; k < nKemp; ++k)
    qZj.col(Kemp(k)).setZero();

  return -logZzj.sum();
}

GTEST_TEST(TestConvexHull, testRandomConvexCombinations) {
  // Generate a set of points, then find a random convex combination of those
  // points, and verify that it's correctly reported as being inside the
  // convex hull
  for (int i = 2; i < 50; ++i) {
    for (int j = 0; j < 500; ++j) {
      MatrixXd pts = MatrixXd::Random(2, i);
      VectorXd weights = VectorXd::Random(i);
      if (weights.minCoeff() < 0) {
        weights = weights.array() -
                  weights.minCoeff();  // make sure they're all nonnegative
      }
      weights = weights.array() / weights.sum();
      Vector2d q = pts * weights;
      EXPECT_TRUE(inConvexHull(pts, q, 1e-8));
    }
  }
}

void  UpdaterMean::costsToWeights(const VectorXd& costs, string weighting_method, double eliteness, VectorXd& weights) const
{
  weights.resize(costs.size());
  if (weighting_method.compare("PI-BB")==0)
  {
    // PI^2 style weighting: continuous, cost exponention
    double h = eliteness; // In PI^2, eliteness parameter is known as "h"
    double range = costs.maxCoeff()-costs.minCoeff();
    if (range==0)
      weights.fill(1);
    else
      weights = (-h*(costs.array()-costs.minCoeff())/range).exp();
  } 
  else if (weighting_method.compare("CMA-ES")==0 || weighting_method.compare("CEM")==0 )
  {
    // CMA-ES and CEM are rank-based, so we must first sort the costs, and the assign a weight to 
    // each rank.
    VectorXd costs_sorted = costs; 
    std::sort(costs_sorted.data(), costs_sorted.data()+costs_sorted.size());
    // In Python this is more elegant because we have argsort.
    // indices = np.argsort(costs)
    // It is possible to do this with fancy lambda functions or std::pair in C++ too, but  I don't
    // mind writing two for loops instead ;-)
    
    weights.fill(0.0);
    int mu = eliteness; // In CMA-ES, eliteness parameter is known as "mu"
    assert(mu<costs.size());
    for (int ii=0; ii<mu; ii++)
    {
      double cur_cost = costs_sorted[ii];
      for (int jj=0; jj<costs.size(); jj++)
      {
        if (costs[jj] == cur_cost)
        {
          if (weighting_method.compare("CEM")==0)
            weights[jj] = 1.0/mu; // CEM
          else
            weights[jj] = log(mu+0.5) - log(ii+1); // CMA-ES
          break;
        }
      }
      
    }
    // For debugging
    //MatrixXd print_mat(3,costs.size());
    //print_mat.row(0) = costs_sorted;
    //print_mat.row(1) = costs;
    //print_mat.row(2) = weights;
    //cout << print_mat << endl;
  }
  else
  {
    cout << __FILE__ << ":" << __LINE__ << ":WARNING: Unknown weighting method '" << weighting_method << "'. Calling with PI-BB weighting." << endl; 
    costsToWeights(costs, "PI-BB", eliteness, weights);
    return;
  }
  
  // Relative standard deviation of total costs
  double mean = weights.mean();
  double std = sqrt((weights.array()-mean).pow(2).mean());
  double rel_std = std/mean;
  if (rel_std<1e-10)
  {
      // Special case: all costs are the same
      // Set same weights for all.
      weights.fill(1);
  }

  // Normalize weights
  weights = weights/weights.sum();

}

double ave(VectorXd y) {
    const double n(y.size());    
    const double ySum(y.sum());
    const double yAve(ySum/n);
    return(yAve);
}