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

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

Eigen::MatrixXd localWeighting( const Eigen::MatrixXd &W, bool isFull, bool isSymmetric)
{
	int n = W.rows();
	double Ls = (W.count()-n)/2; //count number of edges of the graph

	Eigen::MatrixXd C = Eigen::MatrixXd::Zero(n,n);
	double Ws = 0.5*W.sum();

	Eigen::VectorXd D = W.colwise().sum();

	if (isFull)
	{
		if (isSymmetric)
		{
			computeC_symmetric_full(const_cast<Eigen::MatrixXd &>(W),C,D,n);
		}
		else
		{
			// this is a trick to ensure vectorizatoin and no cache misses! some tranpositions have to be made though
			const_cast<Eigen::MatrixXd &>(W).transposeInPlace();
			computeC_symmetric_full(const_cast<Eigen::MatrixXd &>(W),C,D,n);
			const_cast<Eigen::MatrixXd &>(W).transposeInPlace();
			C.transposeInPlace();
			// the original code use vertical access to rows, but is slower
			//compute_C_asymmetric_full(const_cast<Eigen::MatrixXd &>(W),C,D,n);
		}
	}
	else
	{
		if (isSymmetric)
		{
			computeC_symmetric_sparse(const_cast<Eigen::MatrixXd &>(W),C,D,n);
		}
		else
		{
			compute_C_asymmetric_sparse(const_cast<Eigen::MatrixXd &>(W),C,D,n);
		}
	}

	Eigen::MatrixXd G = ((Ls/Ws)*W).cwiseProduct(C);

	Eigen::VectorXd DG = G.colwise().sum();

	for (int i=0; i<n; i++)
	{
		G.row(i)/=DG(i);
	}
	
	return G;
}

void Plsa::plsa_Train(Eigen::MatrixXd& data,int topics,int Max_Iterations, double tolerance){
    //Variables
    data.transposeInPlace();
    srand((unsigned int) time(0));
    M=data.rows();   // number of visual words
    Ni=data.cols();  // number of images
    K=topics;        // number of Topics
    map<int,double> table_likehood; // likelihood Values
    double offset=1e-7; // small offset to avoid numerical  problems in Pw_z

    //incializacion//
    P_z=Eigen::VectorXd::Ones(K)/K; // P(z);
    Pd_z=Eigen::MatrixXd::Random(Ni,K); // P(d|z)
    Pw_z=Eigen::MatrixXd::Random(M,K); //P(w|z)
    Pd_z=Pd_z.cwiseAbs();
    Pw_z=Pw_z.cwiseAbs();

    // normalize
    Normalize(Pd_z);
    Normalize(Pw_z);

    cout <<"P(d|z) "<< Pd_z.colwise().sum()<<endl;
    cout <<"P(w|z) " <<Pw_z.colwise().sum()<<endl;

    //EM algorithm
    for (int i=0;i<Max_Iterations;i++)
    {
        double dli=0;

        E_step();
        M_step(data);

        table_likehood[i]=loglikehood(data);
        Pw_z=(Pw_z.array()+offset).matrix();

        if (i>0)
        {
            dli=table_likehood[i]-table_likehood[i-1];
            if (dli< tolerance)
            {
              cout<<" finalizo en: "<<i<<endl;
              break;
            }
        }
    }
}

/** solves A * X' = B' for X in place */
void chol_solve_t(Eigen::MatrixXd & A, Eigen::MatrixXd & B) {
  if (A.rows() != B.cols()) {throw std::runtime_error("A.rows() must equal B.cols()");}
  B.transposeInPlace();
  chol_solve(A, B);
  B.transposeInPlace();
}