C++ MatrixXd::array方法代码示例

/**
 * Generates an artificial topographyGrid of size numRows x numCols if no
 * topographic data is available.  Results are dumped into topographyGrid.
 * @param topographyGrid A pointer to a zero-initialized Grid of size
 * numRows x numCols.
 * @param numRows The desired number of non-border rows in the resulting matrix.
 * @param numCols The desired number of non-border cols in the resulting matrix.
 */
void simulatetopographyGrid(Grid* topographyGrid, int numRows, int numCols) {
    Eigen::VectorXd refx = refx.LinSpaced(numCols, -2*M_PI, 2*M_PI);
    Eigen::VectorXd refy = refx.LinSpaced(numRows, -2*M_PI, 2*M_PI);
    Eigen::MatrixXd X = refx.replicate(1, numRows);
    X.transposeInPlace();
    Eigen::MatrixXd Y = refy.replicate(1, numCols);

    // Eigen can only deal with two matrices at a time,
    // so split the computation:
    // topographyGrid = sin(X) * sin(Y) * abs(X) * abs(Y) -pi
    Eigen::MatrixXd absXY = X.cwiseAbs().cwiseProduct(Y.cwiseAbs());
    Eigen::MatrixXd sins = X.array().sin().cwiseProduct(Y.array().sin());
    Eigen::MatrixXd temp;
    temp.resize(numRows, numCols);
    temp = absXY.cwiseProduct(sins);

    // All this work to create a matrix of pi...
    Eigen::MatrixXd pi;
    pi.resize(numRows, numCols);
    pi.setConstant(M_PI);
    temp = temp - pi;
    // And the final result.
    topographyGrid->data.block(border, border, numRows, numCols) =
                              temp.block(0, 0, numRows, numCols);
    // Ignore positive values.
    topographyGrid->data =
            topographyGrid->data.unaryExpr(std::ptr_fun(validateDepth));
    topographyGrid->clearNA();
}

TEST(decimate,hemisphere)
{
  // Load a hemisphere centered at the origin. For each original vertex compute
  // its "perfect normal" (i.e., its position treated as unit vectors).
  // Decimate the model and using the birth indices of the output vertices grab
  // their original "perfect normals" and compare them to their current
  // positions treated as unit vectors. If vertices have not moved much, then
  // these should be similar (mostly this is checking if the birth indices are
  // sane).
  Eigen::MatrixXd V,U;
  Eigen::MatrixXi F,G;
  Eigen::VectorXi J,I;
  // Load example mesh: GetParam() will be name of mesh file
  test_common::load_mesh("hemisphere.obj", V, F);
  // Perfect normals from positions
  Eigen::MatrixXd NV = V.rowwise().normalized();
  // Remove half of the faces
  igl::decimate(V,F,F.rows()/2,U,G,J,I);
  // Expect that all normals still point in same direction as original
  Eigen::MatrixXd NU = U.rowwise().normalized();
  Eigen::MatrixXd NVI;
  igl::slice(NV,I,1,NVI);
  ASSERT_EQ(NVI.rows(),NU.rows());
  ASSERT_EQ(NVI.cols(),NU.cols());
  // Dot product
  Eigen::VectorXd D = (NU.array()*NVI.array()).rowwise().sum();
  Eigen::VectorXd O = Eigen::VectorXd::Ones(D.rows());
  // 0.2 chosen to succeed on 256 face hemisphere.obj reduced to 128 faces
  test_common::assert_near(D,O,0.02);
}

/*Predicts where the tracked bounding box will be in the next frame*/
Eigen::Vector4d bb_predict(Eigen::VectorXd const & bb0,
		Eigen::MatrixXd const & pt0, Eigen::MatrixXd const & pt1) {

	Eigen::MatrixXd of(pt1.rows(), pt1.cols());
	of = pt1.array() - pt0.array();
	double dx = median(of.row(0));
	double dy = median(of.row(1));
	unsigned int matCols = pt0.cols(), pos = 0, len = (matCols * (matCols - 1))
			/ 2;
	Eigen::RowVectorXd d1(len);
	Eigen::RowVectorXd d2(len);
	//calculate euclidean distance
	for (unsigned int h = 0; h < matCols - 1; h++)
		for (unsigned int j = h + 1; j < matCols; j++, pos++) {
			d1(pos) = sqrt(((pt0(0, h) - pt0(0, j)) * (pt0(0, h) - pt0(0, j)))
					+ ((pt0(1, h) - pt0(1, j)) * (pt0(1, h) - pt0(1, j))));

			d2(pos) = sqrt(((pt1(0, h) - pt1(0, j)) * (pt1(0, h) - pt1(0, j)))
					+ ((pt1(1, h) - pt1(1, j)) * (pt1(1, h) - pt1(1, j))));
		}

	Eigen::RowVectorXd ds(len);
	//calculate mean value
	ds = d2.array() / d1.array();
	double s = median(ds);
	double s1 = (0.5 * (s - 1) * bb_width(bb0));
	double s2 = (0.5 * (s - 1) * bb_height(bb0));
	Eigen::Vector4d bb1;
	//
	bb1(0) = bb0(0) - s1 + dx;
	bb1(1) = bb0(1) - s2 + dy;
	bb1(2) = bb0(2) + s1 + dx;
	bb1(3) = bb0(3) + s2 + dy;
	return bb1;
}

// Great circle distance, up to a constant of proportionality equal to 2*R
// where R is the earth's radius
Eigen::MatrixXd greatcirc_dist(const Eigen::MatrixXd &X, const Eigen::MatrixXd &Y) {
    int nr = X.rows();
    int nc = Y.rows();
    Eigen::MatrixXd lon1 = X.col(0).replicate(1,nc) * pi_180;
    Eigen::MatrixXd lat1 = X.col(1).replicate(1,nc) * pi_180;
    Eigen::MatrixXd lon2 = Y.col(0).transpose().replicate(nr,1) * pi_180;
    Eigen::MatrixXd lat2 = Y.col(1).transpose().replicate(nr,1) * pi_180;
    Eigen::MatrixXd dlon = 0.5*(lon2 - lon1);
    Eigen::MatrixXd dlat = 0.5*(lat2 - lat1);
    Eigen::MatrixXd a = dlat.array().sin().square().matrix() +
                        (dlon.array().sin().square() * lat1.array().cos() * lat2.array().cos()).matrix();
    Eigen::MatrixXd c = (a.array()<1.0).select(a.array().sqrt(),Eigen::MatrixXd::Ones(nr,nc)).array().asin();
    return (c); // Instead of (2*R*c) where R = 6378137 is the Earth's radius.
}

// Compute the average contour, by calling compute_contour_vals repeatedly
//
// [[Rcpp::export]]
Eigen::MatrixXd rcppstandardize_rates(const Eigen::VectorXd &tiles, const Eigen::VectorXd &rates,
                                      const Eigen::VectorXd &xseed, const Eigen::VectorXd &yseed,
                                      const Eigen::MatrixXd &marks, const Eigen::VectorXd &nmrks,
                                      const std::string &distm) {
    bool use_weighted_mean = true;
    int nxmrks = nmrks(0);
    int nymrks = nmrks(1);
    Eigen::MatrixXd Zvals = Eigen::MatrixXd::Zero(nymrks,nxmrks);
    Eigen::MatrixXd zvals = Eigen::MatrixXd::Zero(nymrks,nxmrks);
    int niters = tiles.size();
    for ( int i = 0, pos = 0 ; i < niters ; i++ ) {
        int now_tiles = (int)tiles(i);
        Eigen::VectorXd now_rates = rates.segment(pos,now_tiles);
        Eigen::VectorXd now_xseed = xseed.segment(pos,now_tiles);
        Eigen::VectorXd now_yseed = yseed.segment(pos,now_tiles);
        Eigen::MatrixXd now_seeds(now_tiles, 2);
        now_seeds << now_xseed,now_yseed;
        if (use_weighted_mean) {
            compute_contour_vals(zvals,marks,now_rates,now_seeds,distm);
            zvals = zvals.array() - zvals.mean();
        } else {
            now_rates = now_rates.array() - now_rates.mean();
            compute_contour_vals(zvals,marks,now_rates,now_seeds,distm);
        }
        Zvals += zvals;
        pos += now_tiles;
    }
    // Do not divide by niters here but in 'average.eems.contours' instead
    // Zvals = Zvals.array() / niters;
    return Zvals.transpose();
}

Eigen::MatrixXd RtHPIS::magnetic_dipole(Eigen::MatrixXd pos, Eigen::MatrixXd pnt, Eigen::MatrixXd ori) {

    double u0 = 1e-7;
    int nchan;
    Eigen::MatrixXd r, r2, r5, x, y, z, mx, my, mz, Tx, Ty, Tz, lf, temp;

    nchan = pnt.rows();

    // Shift the magnetometers so that the dipole is in the origin
    pnt.array().col(0) -= pos(0);pnt.array().col(1) -= pos(1);pnt.array().col(2) -= pos(2);

    r = pnt.array().square().rowwise().sum().sqrt();

   r2 = r5 = x = y = z = mx = my = mz = Tx = Ty = Tz = lf = Eigen::MatrixXd::Zero(nchan,3);

    for(int i = 0;i < nchan;i++) {
            r2.row(i).array().fill(pow(r(i),2));
            r5.row(i).array().fill(pow(r(i),5));
    }

    for(int i = 0;i < nchan;i++) {
        x.row(i).array().fill(pnt(i,0));
        y.row(i).array().fill(pnt(i,1));
        z.row(i).array().fill(pnt(i,2));
    }
    mx.col(0).array().fill(1);my.col(1).array().fill(1);mz.col(2).array().fill(1);

    Tx = 3 * x.cwiseProduct(pnt) - mx.cwiseProduct(r2);
    Ty = 3 * y.cwiseProduct(pnt) - my.cwiseProduct(r2);
    Tz = 3 * z.cwiseProduct(pnt) - mz.cwiseProduct(r2);

    for(int i = 0;i < nchan;i++) {
        lf(i,0) = Tx.row(i).dot(ori.row(i));
        lf(i,1) = Ty.row(i).dot(ori.row(i));
        lf(i,2) = Tz.row(i).dot(ori.row(i));
    }

    for(int i = 0;i < nchan;i++) {
        for(int j = 0;j < 3;j++) {
            lf(i,j) = u0 * lf(i,j)/(4 * M_PI * r5(i,j));
        }
    }
    return lf;
}

int FeatureTransformationEstimator::consensus3D(Eigen::MatrixXd P, Eigen::MatrixXd Q, Eigen::Isometry3d T, double thresh, Eigen::Array<bool, 1, Eigen::Dynamic> &consensusSet)
{
    int consensus = 0;

    P = T * P.colwise().homogeneous();
    Eigen::MatrixXd norms = (P - Q).colwise().norm();
    consensusSet = norms.array() < thresh;
    consensus = consensusSet.count();

    return consensus;
}

int main(int argc, char * argv[])
{
  using namespace Eigen;
  using namespace std;
  using namespace igl;
  VectorXd D;
  if(!read_triangle_mesh("../shared/beetle.off",V,F))
  {
    cout<<"failed to load mesh"<<endl;
  }
  twod = V.col(2).minCoeff()==V.col(2).maxCoeff();
  bbd = (V.colwise().maxCoeff()-V.colwise().minCoeff()).norm();
  SparseMatrix<double> L,M;
  cotmatrix(V,F,L);
  L = (-L).eval();
  massmatrix(V,F,MASSMATRIX_TYPE_DEFAULT,M);
  const size_t k = 5;
  if(!eigs(L,M,k+1,EIGS_TYPE_SM,U,D))
  {
    cout<<"failed."<<endl;
  }
  U = ((U.array()-U.minCoeff())/(U.maxCoeff()-U.minCoeff())).eval();

  igl::viewer::Viewer viewer;
  viewer.callback_key_down = [&](igl::viewer::Viewer & viewer,unsigned char key,int)->bool
  {
    switch(key)
    {
      default: 
        return false;
      case ' ':
      {
        U = U.rightCols(k).eval();
        // Rescale eigen vectors for visualization
        VectorXd Z = 
          bbd*0.5*U.col(c);
        Eigen::MatrixXd C;
        igl::parula(U.col(c).eval(),false,C);
        c = (c+1)%U.cols();
        if(twod)
        {
          V.col(2) = Z;
        }
        viewer.data.set_mesh(V,F);
        viewer.data.compute_normals();
        viewer.data.set_colors(C);
        return true;
      }
    }
  };
  viewer.callback_key_down(viewer,' ',0);
  viewer.core.show_lines = false;
  viewer.launch();
}

void scaleData(Eigen::MatrixXd& data, double min, double max)
{
  if(min >= max)
    throw OpenANNException("Scaling failed: max has to be greater than min!");
  const double minData = data.minCoeff();
  const double maxData = data.maxCoeff();
  const double dataRange = maxData - minData;
  const double desiredRange = max - min;
  const double scaling = desiredRange / dataRange;
  data = data.array() * scaling + (min - minData * scaling);
}

      normal_fullrank operator/=(const normal_fullrank& rhs) {
        static const char* function =
          "stan::variational::normal_fullrank::operator/=";

        stan::math::check_size_match(function,
                             "Dimension of lhs", dimension_,
                             "Dimension of rhs", rhs.dimension());

        mu_.array() /= rhs.mu().array();
        L_chol_.array() /= rhs.L_chol().array();
        return *this;
      }

IGL_INLINE void igl::ViewerData::set_colors(const Eigen::MatrixXd &C)
{
  using namespace std;
  using namespace Eigen;
  // Ambient color should be darker color
  const auto ambient = [](const MatrixXd & C)->MatrixXd
  {
    return 0.1*C;
  };
  // Specular color should be a less saturated and darker color: dampened
  // highlights
  const auto specular = [](const MatrixXd & C)->MatrixXd
  {
    const double grey = 0.3;
    return grey+0.1*(C.array()-grey);
  };
  if (C.rows() == 1)
  {
    for (unsigned i=0;i<V_material_diffuse.rows();++i)
    {
      V_material_diffuse.row(i) = C.row(0);
    }
    V_material_ambient = ambient(V_material_diffuse);
    V_material_specular = specular(V_material_diffuse);

    for (unsigned i=0;i<F_material_diffuse.rows();++i)
    {
      F_material_diffuse.row(i) = C.row(0);
    }
    F_material_ambient = ambient(F_material_diffuse);
    F_material_specular = specular(F_material_diffuse);
  }
  else if (C.rows() == V.rows())
  {
    set_face_based(false);
    V_material_diffuse = C;
    V_material_ambient = ambient(V_material_diffuse);
    V_material_specular = specular(V_material_diffuse);
  }
  else if (C.rows() == F.rows())
  {
    set_face_based(true);
    F_material_diffuse = C;
    F_material_ambient = ambient(F_material_diffuse);
    F_material_specular = specular(F_material_diffuse);
  }
  else
    cerr << "ERROR (set_colors): Please provide a single color, or a color per face or per vertex.";
  dirty |= DIRTY_DIFFUSE;

}

void test_linear_ring()
{
    trace.beginBlock("linear ring");

    const Domain domain(Point(-5,-5), Point(5,5));

    typedef DiscreteExteriorCalculus<1, 2, EigenLinearAlgebraBackend> Calculus;
    Calculus calculus;
    calculus.initKSpace<Domain>(domain);

    for (int kk=-8; kk<10; kk++) calculus.insertSCell( calculus.myKSpace.sCell(Point(-8,kk), kk%2 == 0 ? Calculus::KSpace::POS : Calculus::KSpace::NEG) );
    for (int kk=-8; kk<10; kk++) calculus.insertSCell( calculus.myKSpace.sCell(Point(kk,10), kk%2 == 0 ? Calculus::KSpace::POS : Calculus::KSpace::NEG) );
    for (int kk=10; kk>-8; kk--) calculus.insertSCell( calculus.myKSpace.sCell(Point(10,kk)) );
    for (int kk=10; kk>-8; kk--) calculus.insertSCell( calculus.myKSpace.sCell(Point(kk,-8)) );
    calculus.updateIndexes();

    {
        trace.info() << calculus << endl;
        Board2D board;
        board << domain;
        board << calculus;
        board.saveSVG("ring_structure.svg");
    }

    const Calculus::PrimalDerivative0 d0 = calculus.derivative<0, PRIMAL>();
    display_operator_info("d0", d0);

    const Calculus::PrimalHodge0 h0 = calculus.hodge<0, PRIMAL>();
    display_operator_info("h0", h0);

    const Calculus::DualDerivative0 d0p = calculus.derivative<0, DUAL>();
    display_operator_info("d0p", d0p);

    const Calculus::PrimalHodge1 h1 = calculus.hodge<1, PRIMAL>();
    display_operator_info("h1", h1);

    const Calculus::PrimalIdentity0 laplace = calculus.laplace<PRIMAL>();
    display_operator_info("laplace", laplace);

    const int laplace_size = calculus.kFormLength(0, PRIMAL);
    const Eigen::MatrixXd laplace_dense(laplace.myContainer);

    for (int ii=0; ii<laplace_size; ii++)
        FATAL_ERROR( laplace_dense(ii,ii) == 2 );

    FATAL_ERROR( laplace_dense.array().rowwise().sum().abs().sum() == 0 );
    FATAL_ERROR( laplace_dense.transpose() == laplace_dense );

    trace.endBlock();
}

	void object::test<6>()
	{
		for (int i = 0;i<10;i++)
		{
			int dim = rand()%1000+2;
			int samplenum1 = rand()%1000+100;
			int samplenum2 = rand()%1000+100;

			Eigen::MatrixXd ha = Eigen::MatrixXd::Random(dim,samplenum1);
			Eigen::MatrixXd hb = Eigen::MatrixXd::Random(dim,samplenum2);



			Eigen::VectorXd haa = (ha.array()*ha.array()).colwise().sum();
			Eigen::VectorXd hbb = (hb.array()*hb.array()).colwise().sum();

			Eigen::MatrixXd hdist = -2*ha.transpose()*hb;

			hdist.colwise() += haa;

			hdist.rowwise() += hbb.transpose();

			
			Matrix<double> ga(ha),gb(hb);

			Matrix<double> gdist = -2*ga.transpose()*gb;

			Vector<double> gaa = (ga.array()*ga.array()).colwise().sum();
			Vector<double> gbb = (gb.array()*gb.array()).colwise().sum();

			gdist.colwise() += gaa;
			gdist.rowwise() += gbb;

			ensure(check_diff(hdist,gdist));
		}
	}

dipError RtHPIS::dipfitError(Eigen::MatrixXd pos, Eigen::MatrixXd data, struct sens sensors)
{
    // Variable Declaration
    struct dipError e;
    Eigen::MatrixXd lf, dif;

    // Compute lead field for a magnetic dipole in infinite vacuum
    lf = ft_compute_leadfield(pos, sensors);

    e.moment = pinv(lf) * data;

    dif = data - lf * e.moment;

    e.error = dif.array().square().sum()/data.array().square().sum();

    return e;
}

void KMeansTestCase::clustering()
{
    const int N = 1000;
    const int D = 10;
    Eigen::MatrixXd X(N, D);
    OpenANN::RandomNumberGenerator rng;
    rng.fillNormalDistribution(X);

    OpenANN::KMeans kmeans(D, 5);
    const int batchSize = 200;
    double averageDistToCenter = std::numeric_limits<double>::max();
    for(int i = 0; i < N/batchSize; i++)
    {
        // Data points will be closer to centers after each update
        Eigen::MatrixXd Y = kmeans.fitPartial(X.block(i*batchSize, 0, batchSize, D)).transform(X);
        const double newDistance = Y.array().rowwise().maxCoeff().sum();
        ASSERT(newDistance < averageDistToCenter);
        averageDistToCenter = newDistance;
    }
}