C++ MatrixXi::row方法代码示例

void NeighbourJoining::calcNewD(MatrixXf& currentD, MatrixXi& rowsID, const Pair& p) {
	//calculates distances to new node
	int j = 0;
	for (int i = 0; i < numCurrentNodes - 1; ++i) {
		if (i == p.i)
			j++;
		currentD(numCurrentNodes, i) = (currentD(p.i, i + j) + currentD(p.j, i + j) - currentD(p.i, p.j)) / 2;
		currentD(i, numCurrentNodes) = currentD(numCurrentNodes, i);
	}
	//cout << "distances to new node: " << currentD.row(numCurrentNodes).head(numCurrentNodes-1) <<endl;

	//swaps rows and columns so that the closest pair nodes go right and at the bottom of the matrix
	currentD.row(p.i).head(numCurrentNodes - 1).swap(
			currentD.row(numCurrentNodes - 1).head(numCurrentNodes - 1));

	currentD.col(p.i).head(numCurrentNodes - 1).swap(
			currentD.col(numCurrentNodes - 1).head(numCurrentNodes - 1));

	currentD.row(p.j).head(numCurrentNodes - 1).swap(
			currentD.row(numCurrentNodes).head(numCurrentNodes - 1));

	currentD.col(p.j).head(numCurrentNodes - 1).swap(
			currentD.col(numCurrentNodes).head(numCurrentNodes - 1));

	currentD.diagonal().setZero();
	//cout << "new Matrix:" << endl;  	printMatrix(currentD);

	//adjusts node IDs to new matrix indices
	int newNode = 2 * numObservableNodes - numCurrentNodes;
	rowsID.row(p.i).swap(rowsID.row(numCurrentNodes - 1));
	rowsID.row(p.j).swap(rowsID.row(newNode));

	//cout << "rowsID:   " << rowsID.transpose(); cout << endl;
}

void set_submatrix(MatrixXi& M, const MatrixXi& Msub, const VectorXi& ind)
{
	int k = 0;
	for(int i=0;i<ind.size();i++){
		if(ind(i) == 1){
			M.row(i) = Msub.row(k);
			k++;
		}
	}
}

void WeightedSampling(int m, int N, const MatrixXi &resIndex, vector<int> &sample, int h) {
  sample.reserve(m);
  vector<int> seedIndex;
  RandomSampling(1, N, seedIndex);
  sample[0] = seedIndex[0];
  int prevSelected = seedIndex[0];
  ArrayXf w(N);
  w.setConstant(1);
  for (int i = 1; i < m; i++) {
    ArrayXf new_w(N);
    computeIntersection(resIndex.row(prevSelected), resIndex, h, new_w);
    w(prevSelected) = 0;
    w *= new_w;
    if (w.sum() > 0) {
      map<double, int> cumulative;
      typedef std::map<double, int>::iterator it;
      double acc = 0;
      for (int j = 0; j < N; j++) {
        acc += w(j);
        cumulative[acc] = j; 
      }
      double linear = rand()*acc/RAND_MAX;
      prevSelected = cumulative.upper_bound(linear)->second;
      sample[i] = prevSelected;
    } else {
      vector<int> temp_sample;
      RandomSampling(1, N, temp_sample);
      prevSelected = temp_sample[0];
      sample[i] = prevSelected;
    }
  }
}

void Utils::saveImage(MatrixXi rawData, string outputFile)
{
    QImage* img = new QImage(rawData.cols(), rawData.rows(), QImage::Format_RGB16);
    for (int y = 0; y < img->height(); y++)
    {
        VectorXi a = rawData.row(y);
        memcpy(img->scanLine(y), (void*)&a, img->bytesPerLine());
    }
    QString file = QString::fromUtf8(outputFile.c_str());
    img->save(file);
}

//TODO: Clean up function params size Fixed and size Move are not needed
void Simulation::reIndexTVandTT(
	vector<int> newVertsIndices, 
	int sizeFixed, 
	int sizeMove, 
	MatrixXd& TV, 
	MatrixXi& TT, 
	VectorXd& force, 
	MatrixXd& newTV, 
	MatrixXi& newTT, 
	VectorXd& new_force){
	//apply re-index to TV
	for(unsigned int i=0; i<newVertsIndices.size(); i++){
		newTV.row(i) = TV.row(newVertsIndices[i]);
		new_force.segment<3>(3*i) = force.segment<3>(3*newVertsIndices[i]);
	}
	//create map out of newVertsIndex
	//map keys = newVertsIndex values = old indices in TV
	//map vals = newVertsIndex index = new indices in TV
	map<int, int> oldToNewmap;
	pair<map<int, int>::iterator, bool> err;
	for(unsigned int i=0; i<newVertsIndices.size(); i++){
		err = oldToNewmap.insert(pair<int, int>(newVertsIndices[i], i));
		if(err.second==false){
			cout<<"ERROR::Simulation.cpp::reIndexTVandTT::>>Map already contains this value(";
			cout<< err.first->second <<". Indices should not be repeated"<<endl;
		}
	}

	//Check map, see if its working
	// map<int,int>::iterator it = oldToNewmap.begin();
	// for (it=oldToNewmap.begin(); it!=oldToNewmap.end(); ++it)
	// 	cout << it->first << " => " << it->second << '\n';


	//apply re-index to TT
	for(int i=0; i< TT.rows(); i++){
		for(int j=0; j<4; j++){
			newTT.row(i)[j] = oldToNewmap.find(TT.row(i)[j])->second;
		}
	}
}

MatrixXi get_submatrix(const MatrixXi& M, const VectorXi& ind)
{
	MatrixXi Msub(ind.sum(),M.cols());
	int k = 0;
	for(int i=0;i<M.rows();i++){
		if(ind(i) == 1){
			Msub.row(k) = M.row(i);
			k++;
		}
	}
	return Msub;
}

parameters::parameters(datafile dat, model nv_mod, model ref_mod,parameters ref_param,int compo,int iter){
	const MatrixXi & omega=nv_mod.Get_model(),ref_omega=ref_mod.Get_model(),mat=dat.Get_mat_datafile();
	const int g=omega.rows(),unique=mat.rows();
	m_proba=ref_param.m_proba;
	m_proba_block.resize(g);
	m_param.resize(g);
	for (int k=0;k<g;k++){
		if (k!=compo){
			m_param[k].resize(omega.rowwise().maxCoeff()(k)+1);
			m_param[k]=ref_param.m_param[k];
			m_proba_block[k].resize(unique,omega.rowwise().maxCoeff()(k)+1);
			m_proba_block[k]=ref_param.m_proba_block[k];
			for (int b=0;b<(omega.rowwise().maxCoeff()(k)+1) ;b++){
				if ((omega.row(k).array()==b).any()){
					m_param[k][b]=ref_param.m_param[k][b];
				}
			}
		}else{
			m_param[k].resize(omega.rowwise().maxCoeff()(k)+1);
			m_proba_block[k].resize(unique,omega.rowwise().maxCoeff()(k)+1);
			for (int b=0;b<(omega.rowwise().maxCoeff()(k)+1) ;b++){
				if ((omega.row(k).array()==b).any()){
					if ((((omega.row(k).array()==b)==(ref_omega.row(k).array()==b)).prod())==1){
						m_param[k][b]=ref_param.m_param[k][b];
					}else{
						m_param[k][b]=param_block(k,b,dat,nv_mod,m_proba.col(k).array()/m_proba.rowwise().sum().array(),1);
						if ((omega.row(k).array()==b).count()>1){
							int prem=0;
							while(omega(k,prem)!=b){prem++;}
							if (mat.col(prem).maxCoeff()>5){
									m_param[k][b]=m_param[k][b].Optimise_gamma(k,b,dat,nv_mod,5,m_proba.col(k).array()/m_proba.rowwise().sum().array(),dat.Get_eff_datafile());
							}
						}
					}
				}
			}
		}

	}
	m_propor=uniforme(g);
	Probapost( nv_mod , mat );
	Compte_nbparam(dat,nv_mod);
	Likelihood(dat.Get_eff_datafile());
	Estimation(1,0,iter,dat,nv_mod);
}

TEST(TestEigenHelper, TestSelectRows) {
  unsigned int rows = 4;
  unsigned int cols = 5;
  MatrixXi v(rows, cols);
  v << 1, 2, 3, 4, 5,
       6, 7, 8, 9, 10,
       11, 12, 13, 14, 15,
       16, 17, 18, 19, 20;
  Ints myrows = {1,3};
  MatrixXi z = select_rows<int, Ints >(v, myrows);
  EXPECT_EQ(v.row(1), z.row(0));
  EXPECT_EQ(v.row(3), z.row(1));

  MatrixXd q = v.cast<double>() / 2;
  Ints more_rows = {1,0,2};
  MatrixXd t = select_rows<double, Ints >(q, more_rows);
  for(unsigned int i = 0; i < cols; i++) {
    EXPECT_DOUBLE_EQ(q(1, i), t(0, i));
    EXPECT_DOUBLE_EQ(q(0, i), t(1, i));
    EXPECT_DOUBLE_EQ(q(2, i), t(2, i));
  }

}

void computeIntersection(const MatrixXi &thisIndex, const MatrixXi &index, int h, ArrayXf &new_w) {
  int dataSize = index.rows();
  int M = index.cols();
  for (int i = 0 ; i < dataSize; i++) 
    new_w(i) = intersect(thisIndex.transpose(), index.row(i).transpose(), M, h);
}

IGL_INLINE void igl::reorient_facets_raycast(
  const Eigen::PlainObjectBase<DerivedV> & V,
  const Eigen::PlainObjectBase<DerivedF> & F,
  int rays_total,
  int rays_minimum,
  bool facet_wise,
  bool use_parity,
  bool is_verbose,
  Eigen::PlainObjectBase<DerivedI> & I,
  Eigen::PlainObjectBase<DerivedC> & C)
{
  using namespace Eigen;
  using namespace std;
  assert(F.cols() == 3);
  assert(V.cols() == 3);
  
  // number of faces
  const int m = F.rows();
  
  MatrixXi FF = F;
  if (facet_wise) {
    C.resize(m);
    for (int i = 0; i < m; ++i) C(i) = i;
  
  } else {
    if (is_verbose) cout << "extracting patches... ";
    bfs_orient(F,FF,C);
  }
  if (is_verbose) cout << (C.maxCoeff() + 1)  << " components. ";
  
  // number of patches
  const int num_cc = C.maxCoeff()+1;
  
  // Init Embree
  EmbreeIntersector ei;
  ei.init(V.template cast<float>(),FF);
  
  // face normal
  MatrixXd N;
  per_face_normals(V,FF,N);
  
  // face area
  Matrix<typename DerivedV::Scalar,Dynamic,1> A;
  doublearea(V,FF,A);
  double area_total = A.sum();
  
  // determine number of rays per component according to its area
  VectorXd area_per_component;
  area_per_component.setZero(num_cc);
  for (int f = 0; f < m; ++f)
  {
    area_per_component(C(f)) += A(f);
  }
  VectorXi num_rays_per_component(num_cc);
  for (int c = 0; c < num_cc; ++c)
  {
    num_rays_per_component(c) = max<int>(static_cast<int>(rays_total * area_per_component(c) / area_total), rays_minimum);
  }
  rays_total = num_rays_per_component.sum();
  
  // generate all the rays
  if (is_verbose) cout << "generating rays... ";
  uniform_real_distribution<float> rdist;
  mt19937 prng;
  prng.seed(time(nullptr));
  vector<int     > ray_face;
  vector<Vector3f> ray_ori;
  vector<Vector3f> ray_dir;
  ray_face.reserve(rays_total);
  ray_ori .reserve(rays_total);
  ray_dir .reserve(rays_total);
  for (int c = 0; c < num_cc; ++c)
  {
    if (area_per_component[c] == 0)
    {
      continue;
    }
    vector<int> CF;     // set of faces per component
    vector<double> CF_area;
    for (int f = 0; f < m; ++f)
    {
      if (C(f)==c)
      {
        CF.push_back(f);
        CF_area.push_back(A(f));
      }
    }
    // discrete distribution for random selection of faces with probability proportional to their areas
    discrete_distribution<int> ddist(CF.size(), 0, CF.size(), [&](double i){ return CF_area[static_cast<int>(i)]; });       // simple ctor of (Iter, Iter) not provided by the stupid VC11/12
    for (int i = 0; i < num_rays_per_component[c]; ++i)
    {
      int f = CF[ddist(prng)];          // select face with probability proportional to face area
      float s = rdist(prng);            // random barycentric coordinate (reference: Generating Random Points in Triangles [Turk, Graphics Gems I 1990])
      float t = rdist(prng);
      float sqrt_t = sqrtf(t);
      float a = 1 - sqrt_t;
      float b = (1 - s) * sqrt_t;
      float c = s * sqrt_t;
      Vector3f p = a * V.row(FF(f,0)).template cast<float>().eval()       // be careful with the index!!!
                 + b * V.row(FF(f,1)).template cast<float>().eval()
//.........这里部分代码省略.........

void igl::collapse_small_triangles(
  const Eigen::MatrixXd & V,
  const Eigen::MatrixXi & F,
  const double eps,
  Eigen::MatrixXi & FF)
{
  using namespace Eigen;
  using namespace std;

  // Compute bounding box diagonal length
  double bbd = bounding_box_diagonal(V);
  MatrixXd l;
  edge_lengths(V,F,l);
  VectorXd dblA;
  doublearea(l,dblA);

  // Minimum area tolerance
  const double min_dblarea = 2.0*eps*bbd*bbd;

  Eigen::VectorXi FIM = colon<int>(0,V.rows()-1);
  int num_edge_collapses = 0;
  // Loop over triangles
  for(int f = 0;f<F.rows();f++)
  {
    if(dblA(f) < min_dblarea)
    {
      double minl = 0;
      int minli = -1;
      // Find shortest edge
      for(int e = 0;e<3;e++)
      {
        if(minli==-1 || l(f,e)<minl)
        {
          minli = e;
          minl = l(f,e);
        }
      }
      double maxl = 0;
      int maxli = -1;
      // Find longest edge
      for(int e = 0;e<3;e++)
      {
        if(maxli==-1 || l(f,e)>maxl)
        {
          maxli = e;
          maxl = l(f,e);
        }
      }
      // Be sure that min and max aren't the same
      maxli = (minli==maxli?(minli+1)%3:maxli);

      // Collapse min edge maintaining max edge: i-->j
      // Q: Why this direction?
      int i = maxli;
      int j = ((minli+1)%3 == maxli ? (minli+2)%3: (minli+1)%3);
      assert(i != minli);
      assert(j != minli);
      assert(i != j);
      FIM(F(f,i)) = FIM(F(f,j));
      num_edge_collapses++;
    }
  }

  // Reindex faces
  MatrixXi rF = F;
  // Loop over triangles
  for(int f = 0;f<rF.rows();f++)
  {
    for(int i = 0;i<rF.cols();i++)
    {
      rF(f,i) = FIM(rF(f,i));
    }
  }

  FF.resize(rF.rows(),rF.cols());
  int num_face_collapses=0;
  // Only keep uncollapsed faces
  {
    int ff = 0;
    // Loop over triangles
    for(int f = 0;f<rF.rows();f++)
    {
      bool collapsed = false;
      // Check if any indices are the same
      for(int i = 0;i<rF.cols();i++)
      {
        for(int j = i+1;j<rF.cols();j++)
        {
          if(rF(f,i)==rF(f,j))
          {
            collapsed = true;
            num_face_collapses++;
            break;
          }
        }
      }
      if(!collapsed)
      {
        FF.row(ff++) = rF.row(f);
      }
//.........这里部分代码省略.........

int main(int argc, char *argv[]) {
    cout << "Usage: " << argv[0] << " [FILENAME].[off|obj|ply] [1-7] [sl]"
         << endl;
    cout << "where 1-7 is the cost function to use" << endl;
    cout << "      s = save images at all decimation steps" << endl;
    cout << "      l = disable lighting" << endl;
    cout << endl;
    cout << "Keybindings:" << endl;
    cout << "  [space]  toggle animation." << endl;
    cout << "  'r'  reset." << endl;
    cout << "  '1'  edge collapse." << endl;
    cout << "  '2'  vertex split." << endl;
    cout << "  's'  save screenshot." << endl;
    cout << "  'c'  switch color mode." << endl;
    cout << "  'f'  cycle cost function." << endl;
    cout << endl;
    // Load a closed manifold mesh
    string filename;
    if (argc >= 2) {
        filename = argv[1];
    } else {
        return 0;
    }
    if (argc >= 3) {
        int idx = stoi(argv[2]) - 1;
        cost_function_n = idx;
        if (idx >= 0 && idx < cost_functions.size())
            shortest_edge_and_midpoint = *(cost_functions.begin() + idx);
    }

    if (!igl::read_triangle_mesh(filename, OV, OF)) {
        cout << "could not read mesh from \"" << filename << "\"" << endl;
        return 1;
    }

    // compute normals
    igl::per_face_normals(OV, OF, normals);

    // Prepare array-based edge data structures and priority queue
    // EMAP is a map from faces to edges.
    // Index into it like EMAP(face + i*F.rows()) where i is an edge index
    // between 0 and 2 corresponding to the three edges of a triangle.
    VectorXi EMAP;

    // E is a map from edges to vertices. Given some edge index e,
    // E(e, 0) and E(e, 1) are the two vertices that the edge is composed of.
    MatrixXi E;

    // EF is a map from edges to faces. For some edge index e,
    // EF(e, 0) and E(e, 1) are the two faces that contain the edge e.
    MatrixXi EF;

    // EI is a map from edges to face corners. For some edge index e,
    // EI(e, 0) is the index i such that EMAP(EF(e, 0) + i*F.rows()) == e and
    // EI(e, 1) is the index i such that EMAP(EF(e, 1) + i*F.rows()) == e.
    MatrixXi EI;

    typedef std::set<std::pair<double, int>> PriorityQueue;
    // Q stores the list of possible edge collapses and their costs
    PriorityQueue Q;
    std::vector<PriorityQueue::iterator> Qit;
    // If an edge were collapsed, we'd collapse it to these points:
    MatrixXd C;

    // Keep some info on edge collapses for reversal and debug reasons
    int num_collapsed;
    std::vector<MeshModification> mods;
    std::vector<int> iters;
    int total_decimations = 0;

    const auto &reset_view = [&]() {
        viewer.data.clear();
        viewer.data.set_mesh(V, F);
        switch (color_mode) {
        case DISTANCE_VISUALIZATION:
            generate_distance_field();
        case COST_VISUALIZATION:
            viewer.data.set_colors(colors);
            break;
        case SOLID:
            viewer.data.set_colors(RowVector3d(1.0, 1.0, 1.0));
            break;
        }
        viewer.data.set_face_based(false);
    };

    // Function to reset original mesh and data structures
    const auto &reset = [&]() {
        total_decimations = 0;
        mods.clear();
        iters.clear();
        F = OF;
        V = OV;
        igl::edge_flaps(F, E, EMAP, EF, EI);
        Qit.resize(E.rows());

        C.resize(E.rows(), V.cols());
        colors.resize(V.rows(), 3);
        colors.setZero();
        VectorXd costs(V.rows());
//.........这里部分代码省略.........

RcppExport SEXP BMEclustering(SEXP mat, SEXP moda, SEXP nb_cluster, SEXP partition_initiale, SEXP nb_init,SEXP stop_criterion){
	srand(time(0));

	MatrixXi data=convertMatrix<MatrixXi,NumericMatrix>(mat);
	VectorXi modalite=convertvector<VectorXi,NumericVector>(moda);
	datafile dataF(data,modalite);

	NumericMatrix red=convertMatrix<NumericMatrix,MatrixXi>(dataF.Get_mat_datafile());
	VectorXi partition_vbles=convertvector<VectorXi,NumericVector>(partition_initiale);
	MatrixXi m;
	int g=as<int>(nb_cluster);
	//int borne=as<int>(nbiter);
	m.resize(g,partition_vbles.rows());
	for (int k=0;k<g;k++){
		m.row(k)=partition_vbles;
	}
	NumericMatrix test=convertMatrix<NumericMatrix,MatrixXi>(m);
	int nbinit=as<int>(nb_init);
	int borne=as<int>(stop_criterion);
	MCMCAlgo ref(dataF,m,m.rows(),2,1,2,6,borne);
	for (int ini=0;ini<nbinit;ini++){
		MCMCAlgo test(dataF,m,m.rows(),2,1,2,6,borne);
		if (test.Get_best_bic()>ref.Get_best_bic()){
			ref=test;
		}
	}
	//sauvegarde des caractéristiques du modèle
	NumericMatrix model=convertMatrix<NumericMatrix,MatrixXi>(ref.Get_best_omega());
	double bic=ref.Get_best_bic();
	double likelihood=ref.Get_best_like();
	NumericMatrix probapost=convertMatrix<NumericMatrix,MatrixXd>(ref.Sauv_probapost());
	NumericVector localise=convertvector<NumericVector,VectorXi>(dataF.Get_localise());

	vector< vector< NumericVector > > tau;
	vector< vector< vector< NumericVector > > > delta;
	vector< vector< vector< NumericVector > > > alpha;
	vector< vector< double > > rho;
	tau.resize(g);
	rho.resize(g);
	delta.resize(g);
	alpha.resize(g);
	for (int k=0;k<g;k++){
		tau[k].resize(ref.Get_best_omega().row(k).maxCoeff()+1);
		rho[k].resize(ref.Get_best_omega().row(k).maxCoeff()+1);
		delta[k].resize(ref.Get_best_omega().row(k).maxCoeff()+1);
		alpha[k].resize(ref.Get_best_omega().row(k).maxCoeff()+1);
		for (int b=0;b<=ref.Get_best_omega().row(k).maxCoeff();b++){
		//for (int b=0;b<2;b++){
			//vectblock[k][b]=ref.Get_rho(k,b);
			//vector < VectorXd > passe=ref.Get_alpha(k,b);
			if ((ref.Get_best_omega().row(k).array()==b).count()>0){
				vector<VectorXd> passe=ref.Get_alpha(k,b);
				alpha[k][b].resize((ref.Get_best_omega().row(k).array()==b).count());
				for (int loc=0;loc<((ref.Get_best_omega().row(k).array()==b).count());loc++){
					alpha[k][b][loc]=convertvector<NumericVector,VectorXd>(passe[loc]);
				}
			}
			if ((ref.Get_best_omega().row(k).array()==b).count()>1){
				MatrixXi passe=ref.Get_delta(k,b);
				delta[k][b].resize(passe.cols());
				tau[k][b]=convertvector<NumericVector,VectorXd>(ref.Get_tau(k,b));
				for (int loc=0;loc<passe.cols();loc++){
					delta[k][b][loc]=convertvector<NumericVector,VectorXi>(passe.col(loc));
				}
				//delta[k][b]=convertMatrix<NumericMatrix,MatrixXi>(ref.Get_delta(k,b));
				rho[k][b]=ref.Get_rho(k,b);
			}
		}
	}

	List param=List::create(Rcpp::Named("tau")=tau,Rcpp::Named("rho")=rho,Rcpp::Named("delta")=delta,Rcpp::Named("alpha")=alpha,Rcpp::Named("proportions")=convertvector<NumericVector,VectorXd>(ref.Get_propor()));
    List desc_model = List::create(Rcpp::Named("sigma")=model,Rcpp::Named("bic")=bic,Rcpp::Named("likelihood")=likelihood,Rcpp::Named("probapost")=probapost,Rcpp::Named("partition")=localise,Rcpp::Named("nbcluster")=nb_cluster,Rcpp::Named("parameters")=param);

    return desc_model;
}