C++ Matrix::resize方法代码示例

TEST(AgradRevErrorHandlingMatrix, checkNonzeroSizeMatrix) {
  using stan::agrad::var;
  Eigen::Matrix<var,Eigen::Dynamic,Eigen::Dynamic> y;
  using stan::math::check_nonzero_size;
  var result;
  
  y.resize(3,3);
  EXPECT_TRUE(check_nonzero_size("checkNonzeroSize", "y", y));
  y.resize(2, 3);
  EXPECT_TRUE(check_nonzero_size("checkNonzeroSize", "y", y));

  y.resize(0,0);
  EXPECT_THROW_MSG(check_nonzero_size("checkNonzeroSize", "y", y),
                   std::invalid_argument, "y has size 0");

  std::vector<var> a;
  a.push_back(3.0);
  a.push_back(3.0);
  a.push_back(3.0);
  a.push_back(3.0);


  EXPECT_TRUE(stan::math::check_nonzero_size("checkNonzeroSize",
                                                       "a", a));

  a.resize(2);
  EXPECT_TRUE(stan::math::check_nonzero_size("checkNonzeroSize",
                                                       "a", a));

  a.resize(0);
  EXPECT_THROW_MSG(stan::math::check_nonzero_size("checkNonzeroSize", "a", a),
                   std::invalid_argument,
                   "a has size 0");
}

TEST(MathErrorHandlingMatrix, checkNonzeroSizeMatrix) {
  Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> y;
  double result;
  
  y.resize(3,3);
  EXPECT_TRUE(stan::math::check_nonzero_size("checkNonzeroSize(%1%)",
                                             y, "y", &result));
  y.resize(2, 3);
  EXPECT_TRUE(stan::math::check_nonzero_size("checkNonzeroSize(%1%)",
                                             y, "y", &result));

  y.resize(0,0);
  EXPECT_THROW(stan::math::check_nonzero_size("checkNonzeroSize(%1%)",y, "y", 
                                              &result), 
               std::domain_error);

  std::vector<double> a;
  a.push_back(3.0);
  a.push_back(3.0);
  a.push_back(3.0);
  a.push_back(3.0);


  EXPECT_TRUE(stan::math::check_nonzero_size("checkNonzeroSize(%1%)",
                                             a, "a", &result));

  a.resize(2);
  EXPECT_TRUE(stan::math::check_nonzero_size("checkNonzeroSize(%1%)",
                                             a, "a", &result));

  a.resize(0);
  EXPECT_THROW(stan::math::check_nonzero_size("checkNonzeroSize(%1%)",a, "a", 
                                              &result), 
               std::domain_error);
}

 void
 jacobian(const F& f,
          const Eigen::Matrix<double, Dynamic, 1>& x,
          Eigen::Matrix<double, Dynamic, 1>& fx,
          Eigen::Matrix<double, Dynamic, Dynamic>& J) {
   using Eigen::Matrix;
   using stan::agrad::var;
   start_nested();
   try {
     Matrix<var, Dynamic, 1> x_var(x.size());
     for (int k = 0; k < x.size(); ++k)
       x_var(k) = x(k);
     Matrix<var, Dynamic, 1> fx_var = f(x_var);
     fx.resize(fx_var.size());
     for (int i = 0; i < fx_var.size(); ++i)
       fx(i) = fx_var(i).val();
     J.resize(x.size(), fx_var.size());
     for (int i = 0; i < fx_var.size(); ++i) {
       if (i > 0)
         set_zero_all_adjoints();
       grad(fx_var(i).vi_);
       for (int k = 0; k < x.size(); ++k)
         J(k, i) = x_var(k).adj();
     }
   } catch (const std::exception& e) {
     stan::agrad::recover_memory_nested();
     throw;
   }
   stan::agrad::recover_memory_nested();
 }

int main() {

	gen.seed(static_cast<unsigned int>(std::time(0)));

	const int beam_size = 3;
	const int vocab_size = 10;
	const int start_symbol = 0;
	const int end_symbol = 1;
	const int max_decoding_length = 30;
	const int min_decoding_length = 10;
	const int source_length = 10;
	const precision penalty = 3;

	decoder<precision> d(beam_size,vocab_size,start_symbol,end_symbol,
		max_decoding_length,min_decoding_length,penalty,"");

	Eigen::Matrix<precision,Eigen::Dynamic, Eigen::Dynamic> outputDist;
	Eigen::Matrix<precision, 1, Eigen::Dynamic> normalization;
	outputDist.resize(vocab_size,beam_size);
	normalization.resize(1,beam_size);

	init_output_dist(outputDist,normalization);

	d.init_decoder();

	for(int i=0; i<20; i++) {
		init_output_dist(outputDist,normalization);
		d.expand_hypothesis(outputDist);
		d.print_current_hypotheses();
	}
	d.finish_current_hypotheses(outputDist);
	d.print_current_hypotheses();
}

 void
 hessian(const F& f,
         const Eigen::Matrix<double, Eigen::Dynamic, 1>& x,
         double& fx,
         Eigen::Matrix<double, Eigen::Dynamic, 1>& grad,
         Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>& H) {
   H.resize(x.size(), x.size());
   grad.resize(x.size());
   try {
     for (int i = 0; i < x.size(); ++i) {
       start_nested();
       Eigen::Matrix<fvar<var>, Eigen::Dynamic, 1> x_fvar(x.size());
       for (int j = 0; j < x.size(); ++j)
         x_fvar(j) = fvar<var>(x(j), i == j);
       fvar<var> fx_fvar = f(x_fvar);
       grad(i) = fx_fvar.d_.val();
       if (i == 0) fx = fx_fvar.val_.val();
       stan::math::grad(fx_fvar.d_.vi_);
       for (int j = 0; j < x.size(); ++j)
         H(i, j) = x_fvar(j).val_.adj();
       stan::math::recover_memory_nested();
     }
   } catch (const std::exception& e) {
     stan::math::recover_memory_nested();
     throw;
   }
 }

IGL_INLINE bool igl::png::texture_from_png(
  const std::string png_file,
  Eigen::Matrix<char,Eigen::Dynamic,Eigen::Dynamic>& R,
  Eigen::Matrix<char,Eigen::Dynamic,Eigen::Dynamic>& G,
  Eigen::Matrix<char,Eigen::Dynamic,Eigen::Dynamic>& B,
  Eigen::Matrix<char,Eigen::Dynamic,Eigen::Dynamic>& A
)
{
  int width,height,n;
  unsigned char *data = stbi_load(png_file.c_str(), &width, &height, &n, 4);
  if(data == NULL) {
    return false;
  }

  R.resize(height,width);
  G.resize(height,width);
  B.resize(height,width);
  A.resize(height,width);

  for (unsigned j=0; j<height; ++j) {
    for (unsigned i=0; i<width; ++i) {
      // used to flip with libPNG, but I'm not sure if
      // simply j*width + i wouldn't be better
      // stb_image uses horizontal scanline an starts top-left corner
      R(i,j) = data[4*( (width-1-i) + width * (height-1-j) )];
      G(i,j) = data[4*( (width-1-i) + width * (height-1-j) ) + 1];
      B(i,j) = data[4*( (width-1-i) + width * (height-1-j) ) + 2];
      //A(i,j) = data[4*( (width-1-i) + width * (height-1-j) ) + 3];
    }
  }

  stbi_image_free(data);

  return true;
}

IGL_INLINE void igl::mat_min(
  const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> & X,
  const int dim,
  Eigen::Matrix<T,Eigen::Dynamic,1> & Y,
  Eigen::Matrix<int,Eigen::Dynamic,1> & I)
{
  assert(dim==1||dim==2);

  // output size
  int n = (dim==1?X.cols():X.rows());
  // resize output
  Y.resize(n);
  I.resize(n);

  // loop over dimension opposite of dim
  for(int j = 0;j<n;j++)
  {
    typename Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>::Index PHONY;
    typename Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>::Index i;
    T m;
    if(dim==1)
    {
      m = X.col(j).minCoeff(&i,&PHONY);
    }else
    {
      m = X.row(j).minCoeff(&PHONY,&i);
    }
    Y(j) = m;
    I(j) = i;
  }
}

    Tvector<T> Tmatrix<T>::solve( const Tvector<T>& b, const std::string method ) const
    {
	    assert( ROWS == b.size() ); // Check dimensions are compatible
	    // Convert Tvector to an Eigen matrix
	    Eigen::Matrix<T, -1, 1> B;
	    B.resize(b.size(),1);
	    for (std::size_t i=0; i<b.size(); ++i)
	    {
		    B(i,0) = b.CONTAINER[i];
	    }
	    // Solve the system
	    Eigen::Matrix<T, -1, 1> X;
	    X.resize(b.size(),1);
        assert( method == "LU" || method == "QR" || method == "PartialLU" );
        if ( method == "LU" )
        {
            X = MATRIX.fullPivLu().solve(B);
        }
        if ( method == "PartialLU" )
        {
            X = MATRIX.partialPivLu().solve(B);
        }
        if ( method == "QR" )
        {
            X = MATRIX.fullPivHouseholderQr().solve(B);
        }
	    
	    // Convert back to a Tvector 
	    Tvector<T> x(COLS,0.0);
	    for (std::size_t i=0; i<COLS; ++i)
	    {
		    x.CONTAINER[i] = X(i,0);
	    }	
	    return x;
    }  

int main(){

  ifstream myfile;
  myfile.open("../ex1data1.txt");

  vector<double> vX{};
  vector<double> vy{};

  double d1, d2;
  char c;

  if(myfile.is_open()){
    while(myfile>>d1>>c>>d2){
      vX.push_back(d1);
      vy.push_back(d2);
    }
  }

  int m = vy.size();

  cout<<"Running Gradient Descent ...\n";

  Eigen::Matrix<double,Eigen::Dynamic,2> X;
  Eigen::Matrix<double,Eigen::Dynamic,1> y;

  X.resize(m,2);
  y.resize(m,1);

  for(int i = 0; i < m; ++i){
    X(i,0) = 1;
    X(i,1) = vX[i];
    y(i) = vy[i];
  }

  Eigen::Matrix<double, Eigen::Dynamic, 1> theta = Eigen::MatrixXd::Zero(2,1);

  // Some gradient descent settings
  int iterations = 1500;
  double alpha = 0.01;

  //  gradientDescent(X, y, theta, alpha, iterations);
  gradientDescent(X, y, theta, alpha, iterations);

  cout<<"Theta found by gradient descent: "<<endl;
  cout<<theta<<endl<<endl;

  Eigen::Matrix<double, 1, Eigen::Dynamic> predVec;
  predVec.resize(1,theta.rows());
  predVec<<1, 3.5;
  double predict1 = predVec*theta;

  predVec(1) = 7;
  double predict2 = predVec*theta;

  cout<<"For population = 35,000, we predict a profit of "<<predict1*10000<<endl;
  cout<<"For population = 70,000, we predict a profit of "<<predict2*10000<<endl<<endl;
  return 0;
}

/**
 * Convert from a depth image to 3D camera space coordinates
 * @param depth_image A width x height array of uint16_t depth values
 * @param vertices A width x height array of Vector3f representing the vertices in the depth image in camera space
 * @param normals A width x height array of Vector3f representing the vertices in the depth image in camera space
 */
void Camera::depth_image_to_vertices_and_normals(const uint16_t * depth_image, const uint32_t width, const uint32_t height, Eigen::Matrix<float, 3, Eigen::Dynamic>& vertices, Eigen::Matrix<float, 3, Eigen::Dynamic>& normals ) const {
    using namespace Eigen;

    // Allocate storage for vertx and norms
    vertices.resize( 3, width * height );
    normals.resize( 3, width * height );

    // Run from bottom right corner so we can create normal sin the same pass
    int32_t src_idx = (width * height) - 1;
    for( int16_t y=height-1; y>=0; y-- ) {
        for( int16_t x=width-1; x >= 0; x-- ) {

            // Vertex and normal for this pixel
            Vector3f vertex;
            Vector3f normal { 0, 0, 0 };

            // If this is a valid depth
            uint16_t depth = depth_image[src_idx];
            if( depth != 0 ) {

                // Back project the point into camera 3D space using D(x,y) * Kinv * (x,y,1)T
                Vector2f cam_point = pixel_to_image_plane( x, y );
                vertex = Vector3f { cam_point.x(), cam_point.y(), 1.0f } * depth;

                // Compute normal as v(y,x+1)-v(y,x) cross v(y+1, x)-v(y, x )
                if( (y < static_cast<int16_t>(height - 1 )) && ( x < static_cast<int16_t>(width - 1 ) ) ) {
                    // We have adjacent vertex to right and below that we can extract
                    // Vector[0] is the element to the right of this one
                    // Vector[width] is the element below
                    Vector3f right_neighbour { vertices(0, src_idx + 1),     vertices(1, src_idx + 1),     vertices(2, src_idx + 1) };
                    Vector3f below_neighbour { vertices(0, src_idx + width), vertices(1, src_idx + width), vertices(2, src_idx + width) };

                    // If they are both not BAD
                    if( ( right_neighbour != BAD_VERTEX) && ( below_neighbour != BAD_VERTEX ) ) {
                        right_neighbour -= vertex;
                        below_neighbour -= vertex;

                        // Compute cross product for normal
                        normal = right_neighbour.cross( below_neighbour ).normalized();
                    }
                }
            } else {
                vertex = BAD_VERTEX;
            }

            // Store
            for( int i=0; i<3; i++ ) {
                vertices(i, src_idx) = vertex[i];
                normals(i, src_idx) = normal[i];
            }

            src_idx--;
        }
    }
}

TEST(ErrorHandlingMatrix, checkSquareMatrix) {
  Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> y;
  
  y.resize(3,3);
  EXPECT_TRUE(stan::math::check_square("checkSquareMatrix",
                                                 "y", y));

  y.resize(3, 2);
  EXPECT_THROW(stan::math::check_square("checkSquareMatrix", "y", y), 
               std::invalid_argument);
}

    Tvector<T> Tsparse_matrix<T>::solve( const Tvector<T>& b ) const
    {
        
        assert( ROWS == b.size() );                                 // Check dimensions are compatible
	    // Convert Tvector to an Eigen matrix
	    Eigen::Matrix<T, -1, 1> B;
	    B.resize(b.size(),1);
	    for (std::size_t i=0; i<b.size(); ++i)
	    {
		    B(i,0) = b.CONTAINER[i];
	    }
        
        // Convert Tsparse_matrix to an Eigen::SparseMatrix
        Eigen::SparseMatrix<T> A(ROWS,COLS);  // Declare Eigen sparse matrix
        std::vector<Eigen::Triplet<T>> triplet_list;                // Declare list of triplets 
        for ( std::size_t i=0; i<ROWS; ++i )
        {
            if ( !S_MATRIX[i].isempty() )                           // Check that the row is not empty 
            {   
                std::vector<std::size_t> index;
                std::vector<T> element;
                index = S_MATRIX[i].index_list();
                element = S_MATRIX[i].element_list();
                std::size_t J = index.size();
                for ( std::size_t j=0; j<J; ++j )
                {
                    triplet_list.push_back( Eigen::Triplet<T>( i, index[j], element[j] ));
                }
            }
        } 
        A.setFromTriplets(triplet_list.begin(), triplet_list.end());
        
        // Setup and solve the system
        Eigen::SparseLU< Eigen::SparseMatrix<T> > solverA;
		//Eigen::BiCGSTAB<Eigen::SparseMatrix<T>> solverA;
		//solverA.analyzePattern(A);
		//solverA.factorize(A);
        solverA.compute(A);

        Eigen::Matrix<T, -1, 1> X;
	    X.resize(b.size(),1);
        X = solverA.solve(B);
        
        // Convert back to a Tvector 
	    Tvector<T> x(COLS,0.0);
	    for (std::size_t i=0; i<COLS; ++i)
	    {
		    x.CONTAINER[i] = X(i,0);
	    }	
	    return x;
    }   

TEST(ErrorHandlingMatrix, checkVectorMatrix) {
  Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> x;
  
  x.resize(3,3);
  EXPECT_THROW(stan::math::check_vector("checkVector", "x", x),
               std::invalid_argument);
  x.resize(0,0);
  EXPECT_THROW(stan::math::check_vector("checkVector", "x", x),
               std::invalid_argument);

  x.resize(1,5);
  EXPECT_TRUE(stan::math::check_vector("checkVector", "x", x));

  x.resize(5,1);
  EXPECT_TRUE(stan::math::check_vector("checkVector", "x", x));
}

 static inline void _fast_vector_copy(Eigen::Matrix<T, Eigen::Dynamic, 1>& v_to, const Eigen::Matrix<T, Eigen::Dynamic, 1>& v_from) {
   int sz = v_from.size();
   v_to.resize(sz);
   if (sz > 0) {
     std::memcpy(&v_to(0), &v_from(0), v_from.size() * sizeof(double));
   }
 }

SEXP
vectorToMatVCL(SEXP A, const int nr, const int nc, int device_flag)
{
    Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> Am;
    Am = Rcpp::as<Eigen::Matrix<T, Eigen::Dynamic, 1> >(A);
//    Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> temp = as<Eigen::Matrix<T, Eigen::Dynamic, 1> >(A);
    Am.resize(nr, nc);

//    Eigen::Map<Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> > Amm(&Am(0), nr, nc);

    dynVCLMat<T> *mat = new dynVCLMat<T>(Am, nr, nc, device_flag);
    Rcpp::XPtr<dynVCLMat<T> > pMat(mat);
    return pMat;

//    Eigen::Map<Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> > Amm(&Am(0), nr, nc);
//
//    //use only GPUs:
//    if(device_flag == 0){
//        long id = 0;
//        viennacl::ocl::set_context_device_type(id, viennacl::ocl::gpu_tag());
//    }
//
//    viennacl::matrix<T> *vcl_A = new viennacl::matrix<T>(nr, nc);
//
//    viennacl::copy(Amm, *vcl_A);
//
//    Rcpp::XPtr<viennacl::matrix<T> > pMat(vcl_A);
//    return pMat;
}