C++ SparseMatrix::outerSize方法代码示例

IGL_INLINE void igl::cat(
    const int dim, 
    const Eigen::SparseMatrix<Scalar> & A, 
    const Eigen::SparseMatrix<Scalar> & B, 
    Eigen::SparseMatrix<Scalar> & C)
{
  assert(dim == 1 || dim == 2);
  using namespace Eigen;
  // Special case if B or A is empty
  if(A.size() == 0)
  {
    C = B;
    return;
  }
  if(B.size() == 0)
  {
    C = A;
    return;
  }

  DynamicSparseMatrix<Scalar, RowMajor> dyn_C;
  if(dim == 1)
  {
    assert(A.cols() == B.cols());
    dyn_C.resize(A.rows()+B.rows(),A.cols());
  }else if(dim == 2)
  {
    assert(A.rows() == B.rows());
    dyn_C.resize(A.rows(),A.cols()+B.cols());
  }else
  {
    fprintf(stderr,"cat.h: Error: Unsupported dimension %d\n",dim);
  }

  dyn_C.reserve(A.nonZeros()+B.nonZeros());

  // Iterate over outside of A
  for(int k=0; k<A.outerSize(); ++k)
  {
    // Iterate over inside
    for(typename SparseMatrix<Scalar>::InnerIterator it (A,k); it; ++it)
    {
      dyn_C.coeffRef(it.row(),it.col()) += it.value();
    }
  }

  // Iterate over outside of B
  for(int k=0; k<B.outerSize(); ++k)
  {
    // Iterate over inside
    for(typename SparseMatrix<Scalar>::InnerIterator it (B,k); it; ++it)
    {
      int r = (dim == 1 ? A.rows()+it.row() : it.row());
      int c = (dim == 2 ? A.cols()+it.col() : it.col());
      dyn_C.coeffRef(r,c) += it.value();
    }
  }

  C = SparseMatrix<Scalar>(dyn_C);
}

Eigen::VectorXd KronProdSPMat2(
    Eigen::SparseMatrix<double, Eigen::RowMajor> a0,
    Eigen::SparseMatrix<double, Eigen::RowMajor> a1,
    Eigen::VectorXd y) {

    signal(SIGSEGV, handler);   // install our handler
    Eigen::VectorXd retvec;
    retvec.setZero( a0.rows() * a1.rows()  );

    //loop rows a0
    for (int row_idx0=0; row_idx0<a0.outerSize(); ++row_idx0) {
        int row_offset1 = row_idx0;
        row_offset1    *= a1.rows();

        // loop rows a1
        for (int row_idx1=0; row_idx1<a1.outerSize(); ++row_idx1) {

            // loop cols a0 (non-zero elements only)
            for (Eigen::SparseMatrix<double,RowMajor>::InnerIterator it0(a0,row_idx0); it0; ++it0) {
                int col_offset1 = it0.index();
                col_offset1    *= a1.innerSize();
                double factor1 = it0.value();

                for (Eigen::SparseMatrix<double,RowMajor>::InnerIterator it1(a1,row_idx1); it1; ++it1) {
                    retvec( row_offset1 + row_idx1 ) += factor1 * it1.value() * y( col_offset1 + it1.index() );
                }
            }
        }
    }
    return retvec;
}

IGL_INLINE void igl::matlab::prepare_lhs_double(
    const Eigen::SparseMatrix<Vtype> & M,
    mxArray *plhs[])
{
    using namespace std;
    const int m = M.rows();
    const int n = M.cols();
    // THIS WILL NOT WORK FOR ROW-MAJOR
    assert(n==M.outerSize());
    const int nzmax = M.nonZeros();
    plhs[0] = mxCreateSparse(m, n, nzmax, mxREAL);
    mxArray * mx_data = plhs[0];
    // Copy data immediately
    double * pr = mxGetPr(mx_data);
    mwIndex * ir = mxGetIr(mx_data);
    mwIndex * jc = mxGetJc(mx_data);

    // Iterate over outside
    int k = 0;
    for(int j=0; j<M.outerSize(); j++)
    {
        jc[j] = k;
        // Iterate over inside
        for(typename Eigen::SparseMatrix<Vtype>::InnerIterator it (M,j); it; ++it)
        {
            // copy (cast to double)
            pr[k] = it.value();
            ir[k] = it.row();
            k++;
        }
    }
    jc[M.outerSize()] = k;

}

IGL_INLINE void igl::find(
  const Eigen::SparseMatrix<T>& X,
  Eigen::DenseBase<DerivedI> & I,
  Eigen::DenseBase<DerivedJ> & J,
  Eigen::DenseBase<DerivedV> & V)
{
  // Resize outputs to fit nonzeros
  I.derived().resize(X.nonZeros(),1);
  J.derived().resize(X.nonZeros(),1);
  V.derived().resize(X.nonZeros(),1);

  int i = 0;
  // Iterate over outside
  for(int k=0; k<X.outerSize(); ++k)
  {
    // Iterate over inside
    for(typename Eigen::SparseMatrix<T>::InnerIterator it (X,k); it; ++it)
    {
      V(i) = it.value();
      I(i) = it.row();
      J(i) = it.col();
      i++;
    }
  }
}

void ProbitNoise::evalModel(Eigen::SparseMatrix<double> & Ytest, const int n, Eigen::VectorXd & predictions, Eigen::VectorXd & predictions_var, const Eigen::MatrixXd &cols, const Eigen::MatrixXd &rows, double mean_rating) {
  const unsigned N = Ytest.nonZeros();
  Eigen::VectorXd pred(N);
  Eigen::VectorXd test(N);

// #pragma omp parallel for schedule(dynamic,8) reduction(+:se, se_avg) <- dark magic :)
  for (int k = 0; k < Ytest.outerSize(); ++k) {
    int idx = Ytest.outerIndexPtr()[k];
    for (Eigen::SparseMatrix<double>::InnerIterator it(Ytest,k); it; ++it) {
     pred[idx] = nCDF(cols.col(it.col()).dot(rows.col(it.row())));
     test[idx] = it.value();

      // https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Online_algorithm
      double pred_avg;
      if (n == 0) {
        pred_avg = pred[idx];
      } else {
        double delta = pred[idx] - predictions[idx];
        pred_avg = (predictions[idx] + delta / (n + 1));
        predictions_var[idx] += delta * (pred[idx] - pred_avg);
      }
      predictions[idx++] = pred_avg;

   }
  }
  auc_test_onesample = auc(pred,test);
  auc_test = auc(predictions, test);
}

void ApplyConstraints(Eigen::SparseMatrix<float>& K, const std::vector<Constraint>& constraints)
{
	std::vector<int> indicesToConstraint;

	for (std::vector<Constraint>::const_iterator it = constraints.begin(); it != constraints.end(); ++it)
	{
		if (it->type & Constraint::UX)
		{
			indicesToConstraint.push_back(2 * it->node + 0);
		}
		if (it->type & Constraint::UY)
		{
			indicesToConstraint.push_back(2 * it->node + 1);
		}
	}

	for (int k = 0; k < K.outerSize(); ++k)
	{
		for (Eigen::SparseMatrix<float>::InnerIterator it(K, k); it; ++it)
		{
			for (std::vector<int>::iterator idit = indicesToConstraint.begin(); idit != indicesToConstraint.end(); ++idit)
			{
				SetConstraints(it, *idit);
			}
		}
	}
}

IGL_INLINE void igl::slice_into(
  const Eigen::SparseMatrix<T>& X,
  const Eigen::Matrix<int,Eigen::Dynamic,1> & R,
  const Eigen::Matrix<int,Eigen::Dynamic,1> & C,
  Eigen::SparseMatrix<T>& Y)
{

  int xm = X.rows();
  int xn = X.cols();
  assert(R.size() == xm);
  assert(C.size() == xn);
#ifndef NDEBUG
  int ym = Y.size();
  int yn = Y.size();
  assert(R.minCoeff() >= 0);
  assert(R.maxCoeff() < ym);
  assert(C.minCoeff() >= 0);
  assert(C.maxCoeff() < yn);
#endif

  // create temporary dynamic sparse matrix
  Eigen::DynamicSparseMatrix<T, Eigen::RowMajor>  dyn_Y(Y);
  // Iterate over outside
  for(int k=0; k<X.outerSize(); ++k)
  {
    // Iterate over inside
    for(typename Eigen::SparseMatrix<T>::InnerIterator it (X,k); it; ++it)
    {
      dyn_Y.coeffRef(R(it.row()),C(it.col())) = it.value();
    }
  }
  Y = Eigen::SparseMatrix<T>(dyn_Y);
}

// inserts the sparse matrix 'ins' into the sparse matrix 'original' in the place given by 'row' and 'col' integers
void insertSparseBlock(const Eigen::SparseMatrix<Scalar>& ins, Eigen::SparseMatrix<Scalar>& original, const unsigned int& row, const unsigned int& col)
{
  for (int k=0; k<ins.outerSize(); ++k)
    for (Eigen::SparseMatrix<Scalar>::InnerIterator iti(ins,k); iti; ++iti)
      original.coeffRef(iti.row() + row, iti.col() + col) = iti.value();

  original.makeCompressed();
}

double maxAbsCoeff(const Eigen::SparseMatrix<T> &mat) {
  double maxval;
  int count = 0;
  for (int k = 0; k < mat.outerSize(); ++k) {
    for (typename Eigen::SparseMatrix<T>::InnerIterator it(mat, k); it; ++it) {
      maxval = count == 0 ? std::abs(it.value()) : std::max(std::abs(it.value()), maxval);
    }
  }
  return maxval;
}

IGL_INLINE void igl::components(
  const Eigen::SparseMatrix<AScalar> & A,
  Eigen::PlainObjectBase<DerivedC> & C,
  Eigen::PlainObjectBase<Derivedcounts> & counts)
{
  using namespace Eigen;
  using namespace std;
  assert(A.rows() == A.cols() && "A should be square.");
  const size_t n = A.rows();
  Array<bool,Dynamic,1> seen = Array<bool,Dynamic,1>::Zero(n,1);
  C.resize(n,1);
  typename DerivedC::Scalar id = 0;
  vector<typename Derivedcounts::Scalar> vcounts;
  // breadth first search
  for(int k=0; k<A.outerSize(); ++k)
  {
    if(seen(k))
    {
      continue;
    }
    queue<int> Q;
    Q.push(k);
    vcounts.push_back(0);
    while(!Q.empty())
    {
      const int f = Q.front();
      Q.pop();
      if(seen(f))
      {
        continue;
      }
      seen(f) = true;
      C(f,0) = id;
      vcounts[id]++;
      // Iterate over inside
      for(typename SparseMatrix<AScalar>::InnerIterator it (A,f); it; ++it)
      {
        const int g = it.index();
        if(!seen(g) && it.value())
        {
          Q.push(g);
        }
      }
    }
    id++;
  }
  assert((size_t) id == vcounts.size());
  const size_t ncc = vcounts.size();
  assert((size_t)C.maxCoeff()+1 == ncc);
  counts.resize(ncc,1);
  for(size_t i = 0;i<ncc;i++)
  {
    counts(i) = vcounts[i];
  }
}

IGL_INLINE void igl::repdiag(
  const Eigen::SparseMatrix<T>& A,
  const int d,
  Eigen::SparseMatrix<T>& B)
{
  using namespace std;
  using namespace Eigen;
  int m = A.rows();
  int n = A.cols();

  vector<Triplet<T> > IJV;
  IJV.reserve(A.nonZeros()*d);
  // Loop outer level
  for (int k=0; k<A.outerSize(); ++k)
  {
    // loop inner level
    for (typename Eigen::SparseMatrix<T>::InnerIterator it(A,k); it; ++it)
    {
      for(int i = 0;i<d;i++)
      {
        IJV.push_back(Triplet<T>(i*m+it.row(),i*n+it.col(),it.value()));
      }
    }
  }
  B.resize(m*d,n*d);
  B.setFromTriplets(IJV.begin(),IJV.end());
  

  // Q: Why is this **Very** slow?

  //int m = A.rows();
  //int n = A.cols();

  //B.resize(m*d,n*d);
  //// Reserve enough space for new non zeros
  //B.reserve(d*A.nonZeros());

  //// loop over reps
  //for(int i=0;i<d;i++)
  //{
  //  // Loop outer level
  //  for (int k=0; k<A.outerSize(); ++k)
  //  {
  //    // loop inner level
  //    for (typename Eigen::SparseMatrix<T>::InnerIterator it(A,k); it; ++it)
  //    {
  //      B.insert(i*m+it.row(),i*n+it.col()) = it.value();
  //    }
  //  }
  //}
  //B.makeCompressed();
}

IGL_INLINE void igl::adjacency_matrix(
  const Eigen::PlainObjectBase<DerivedF> & F, 
  Eigen::SparseMatrix<T>& A)
{
  using namespace std;
  using namespace Eigen;
  typedef typename DerivedF::Scalar Index;

  typedef Triplet<T> IJV;
  vector<IJV > ijv;
  ijv.reserve(F.size()*2);
  // Loop over faces
  for(int i = 0;i<F.rows();i++)
  {
    // Loop over this face
    for(int j = 0;j<F.cols();j++)
    {
      // Get indices of edge: s --> d
      Index s = F(i,j);
      Index d = F(i,(j+1)%F.cols());
      ijv.push_back(IJV(s,d,1));
      ijv.push_back(IJV(d,s,1));
    }
  }

  const Index n = F.maxCoeff()+1;
  A.resize(n,n);
  switch(F.cols())
  {
    case 3:
      A.reserve(6*(F.maxCoeff()+1));
      break;
    case 4:
      A.reserve(26*(F.maxCoeff()+1));
      break;
  }
  A.setFromTriplets(ijv.begin(),ijv.end());

  // Force all non-zeros to be one

  // Iterate over outside
  for(int k=0; k<A.outerSize(); ++k)
  {
    // Iterate over inside
    for(typename Eigen::SparseMatrix<T>::InnerIterator it (A,k); it; ++it)
    {
      assert(it.value() != 0);
      A.coeffRef(it.row(),it.col()) = 1;
    }
  }
}

IGL_INLINE void igl::harwell_boeing(
  const Eigen::SparseMatrix<Scalar> & A,
  int & num_rows,
  std::vector<Scalar> & V,
  std::vector<Index> & R,
  std::vector<Index> & C)
{
  num_rows = A.rows();
  int num_cols = A.cols();
  int nnz = A.nonZeros();
  V.resize(nnz);
  R.resize(nnz);
  C.resize(num_cols+1);

  // Assumes outersize is columns
  assert(A.cols() == A.outerSize());
  int column_pointer = 0;
  int i = 0;
  int k = 0;
  // Iterate over outside
  for(; k<A.outerSize(); ++k)
  {
    C[k] = column_pointer;
    // Iterate over inside
    for(typename Eigen::SparseMatrix<Scalar>::InnerIterator it (A,k); it; ++it)
    {
      V[i] = it.value();
      R[i] = it.row();
      i++;
      // Also increment column pointer
      column_pointer++; 
    }
  }
  // by convention C[num_cols] = nnz
  C[k] = column_pointer;
}

Eigen::SparseMatrix<Type> kronecker(Eigen::SparseMatrix<Type> x,
				    Eigen::SparseMatrix<Type> y){
  typedef Eigen::Triplet<Type> T;
  typedef typename Eigen::SparseMatrix<Type>::InnerIterator Iterator;
  std::vector<T> tripletList;
  int n1=x.rows(),n2=x.cols(),n3=y.rows(),n4=y.cols();
  int i,j,k,l;
  // Loop over nonzeros of x
  for (int cx=0; cx<x.outerSize(); cx++)
    for (Iterator itx(x,cx); itx; ++itx)
      // Loop over nonzeros of y
      for (int cy=0; cy<y.outerSize(); cy++)
  	for (Iterator ity(y,cy); ity; ++ity)
  	  {
  	    i=itx.row();
  	    j=itx.col();
  	    k=ity.row();
  	    l=ity.col();
  	    tripletList.push_back(T(i*n3+k,j*n4+l, itx.value()*ity.value() ));
   	  }
  Eigen::SparseMatrix<Type> mat(n1*n3,n2*n4);
  mat.setFromTriplets(tripletList.begin(), tripletList.end());
  return mat;
}

IGL_INLINE void igl::full(
  const Eigen::SparseMatrix<T> & A,
  Eigen::PlainObjectBase<DerivedB>& B)
{
  assert(false && "Obsolete. Just call B = Matrix(A)");
  using namespace Eigen;
  B = PlainObjectBase<DerivedB >::Zero(A.rows(),A.cols());
  // Iterate over outside
  for(int k=0; k<A.outerSize(); ++k)
  {
    // Iterate over inside
    for(typename SparseMatrix<T>::InnerIterator it (A,k); it; ++it)
    {
      B(it.row(),it.col()) = it.value();
    }
  }
}