C++ MatrixXf::diagonal方法代码示例

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;
}

VectorXf param_sensitivity_widget::computeSensitivity(
    MatrixXf &parameterMatrix, VectorXf &responseVector)
{
    MatrixXf Ctemp = parameterMatrix.transpose()*parameterMatrix;
    MatrixXf C;
    C = Ctemp.inverse();

    VectorXf b = C*parameterMatrix.transpose()*responseVector;

    VectorXf Y_hat = parameterMatrix*b;

    int p = b.rows();

    VectorXf sigma2Vec = responseVector-Y_hat;

    float sigma2 = sigma2Vec.squaredNorm();
    sigma2= sigma2/(parameterMatrix.rows() - p);

    Ctemp = C*sigma2;

    MatrixXf denominator = Ctemp.diagonal();

    // Do element-wise division
    VectorXf t = b;
    for (int i = 0; i < b.rows(); i++)
    {
        t(i) = abs(b(i)/sqrt(denominator(i)));
    }

    return t;
}

void NeighbourJoining::calcQ(const MatrixXf& currentD, MatrixXf& Q) {
	//calculates sum of the rows of the distance matrix
	Matrix<float, latentNodes, 1> sumRows;
	//MatrixXf sumRows (numObservableNodes+1,1) ;
	sumRows.head(numCurrentNodes) = currentD.topLeftCorner(numCurrentNodes,
			numCurrentNodes).colwise().sum();
	//cout << sumRows.head(numCurrentNodes); cout << endl;
	//Q = (n-2) * currentD
	Q.topLeftCorner(numCurrentNodes, numCurrentNodes) = (numCurrentNodes - 2)
			* currentD.topLeftCorner(numCurrentNodes, numCurrentNodes);
	//each row in Q - sumRows and each column in Q - sumRows
	Q.topLeftCorner(numCurrentNodes, numCurrentNodes).colwise() -= sumRows.head(
			numCurrentNodes);
	Q.topLeftCorner(numCurrentNodes, numCurrentNodes).rowwise() -= sumRows.head(
			numCurrentNodes).transpose();

	Q.diagonal().setZero();
	//cout << "Q Matrix:" << endl;		printMatrix(Q);
}

void sparsify_soft_relative_to_row_col_max (
    const MatrixXf & dense,
    MatrixSf & sparse,
    float relthresh,
    bool ignore_diagonal)
{
  ASSERT_LT(0, relthresh);

  LOG("sparsifying " << dense.rows() << " x " << dense.cols()
      << " positive matrix to relative threshold " << relthresh);

  VectorXf row_max;
  VectorXf col_max;

  if (ignore_diagonal) {

    VectorXf diag = dense.diagonal();
    const_cast<MatrixXf &>(dense).diagonal().setZero();

    row_max = dense.rowwise().maxCoeff();
    col_max = dense.colwise().maxCoeff();

    const_cast<MatrixXf &>(dense).diagonal() = diag;

  } else {

    row_max = dense.rowwise().maxCoeff();
    col_max = dense.colwise().maxCoeff();

  }

  const int I = dense.rows();
  const int J = dense.cols();

  sparse.resize(I,J);

  double sum_dense = 0;
  double sum_sparse = 0;

  for (int j = 0; j < J; ++j) {
    for (int i = 0; i < I; ++i) {

      const float dense_ij = dense(i,j);
      sum_dense += dense_ij;

      const float thresh = relthresh * min(row_max(i), col_max(j));
      if (dense_ij > thresh) {

        sparse.insert(i,j) = dense_ij;
        sum_sparse += dense_ij;
      }
    }
  }

  sparse.finalize();

  float density = sparse.nonZeros() / float(I * J);
  float loss = (sum_dense - sum_sparse) / sum_dense;
  LOG("sparsifying to density " << density
      << " loses " << (100 * loss) << "% of mass");
}