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

double compute_max_step_from_singularities(const Eigen::MatrixXd& uv,
                                            const Eigen::MatrixXi& F,
                                            Eigen::MatrixXd& d) {
    double max_step = INFINITY;

    // The if statement is outside the for loops to avoid branching/ease parallelizing
    if (uv.cols() == 2) {
      for (int f = 0; f < F.rows(); f++) {
        double min_positive_root = get_min_pos_root_2D(uv,F,d,f);
        max_step = min(max_step, min_positive_root);
      }
    } else { // volumetric deformation
      for (int f = 0; f < F.rows(); f++) {
        double min_positive_root = get_min_pos_root_3D(uv,F,d,f);
        max_step = min(max_step, min_positive_root);
      }
    }
    return max_step;
 }

  /*
  Train the classifier with the dataset of breastcancer patients from the
  LibSVM Libary
  */
  void TrainSVMClassifier_BreastCancerDataSet_shouldReturnTrue()
  {
    /* Declarating an featurematrixdataset, the first matrix
    of the matrixpair is the trainingmatrix and the second one is the testmatrix.*/
    std::pair<MatrixDoubleType,MatrixDoubleType> matrixDouble;
    matrixDouble = convertCSVToMatrix<double>(GetTestDataFilePath("Classification/FeaturematrixBreastcancer.csv"),';',0.5,true);
    m_TrainingMatrixX = matrixDouble.first;
    m_TestXPredict = matrixDouble.second;

    /* The declaration of the labelmatrixdataset is equivalent to the declaration
    of the featurematrixdataset.*/
    std::pair<MatrixIntType,MatrixIntType> matrixInt;
    matrixInt = convertCSVToMatrix<int>(GetTestDataFilePath("Classification/LabelmatrixBreastcancer.csv"),';',0.5,false);
    m_TrainingLabelMatrixY = matrixInt.first;
    m_TestYPredict = matrixInt.second;

    /* Setting of the SVM-Parameters*/
    classifier = mitk::LibSVMClassifier::New();
    classifier->SetGamma(1/(double)(m_TrainingMatrixX.cols()));
    classifier->SetSvmType(0);
    classifier->SetKernelType(2);

    /* Train the classifier, by giving trainingdataset for the labels and features.
    The result in an colunmvector of the labels.*/
    classifier->Train(m_TrainingMatrixX,m_TrainingLabelMatrixY);
    Eigen::MatrixXi classes = classifier->Predict(m_TestXPredict);

    /* Testing the matching between the calculated colunmvector and the result
    of the SVM */
    unsigned int maxrows = classes.rows();

    bool isYPredictVector = false;
    int count = 0;

    for(int i= 0; i < maxrows; i++){
      if(classes(i,0) == m_TestYPredict(i,0)){
        isYPredictVector = true;
        count++;
      }
    }
    MITK_INFO << 100*count/(double)(maxrows) << "%";
    MITK_TEST_CONDITION(isIntervall<int>(m_TestYPredict,classes,75,100),"Testvalue is in range.");
  }

bool BlockSparseMatrix::ConjugateGradientSolve(const Eigen::MatrixXd& rhs, Eigen::MatrixXd& sol)
{
	return LltSolve(rhs, sol);

	MatrixSolver* solver = new PCGSolver();
	int numRows = rhs.rows();
	int numCols = rhs.cols();
	int num = numRows * numCols;
	double* result = new double[num];
	double* rhsData = new double[num];
	int ithElem = 0;
	for (int i = 0; i < numCols; ++i)
	{
		for (int j = 0; j < numRows; ++j)
		{
			rhsData[ithElem++] = rhs(j, i);
		}
	}

	if (mbCSRDirty)
	{
		mCSREquivalent.ConvertFromBlockSparseMatrix(*this);
		mbCSRDirty = false;
	}

	solver->SetMatrix(&mCSREquivalent);
	solver->SetRHS(rhsData, numCols);
	solver->Solve(result, true);
	sol = Eigen::MatrixXd::Zero(numRows, numCols);
	ithElem = 0;
	for (int i = 0; i < numCols; ++i)
	{
		for (int j = 0; j < numRows; ++j)
		{
			sol(j, i) = result[ithElem++];
		}
	}

	delete[] rhsData;
	delete[] result;
	delete solver;
	return true;
}

// Derived from code by Yohann Solaro ( http://listengine.tuxfamily.org/lists.tuxfamily.org/eigen/2010/01/msg00187.html )
// see : http://en.wikipedia.org/wiki/Moore-Penrose_pseudoinverse#The_general_case_and_the_SVD_method
Eigen::MatrixXd pinv( const Eigen::MatrixXd &b, double rcond )
{
  // TODO: Figure out why it wants fewer rows than columns
  // if ( a.rows()<a.cols() )
  // return false;
  bool flip = false;
  Eigen::MatrixXd a;
  if( a.rows() < a.cols() )
  {
    a = b.transpose();
    flip = true;
  }
  else
    a = b;
  // SVD
  Eigen::JacobiSVD<Eigen::MatrixXd> svdA;
  svdA.compute( a, Eigen::ComputeFullU | Eigen::ComputeThinV );
  Eigen::JacobiSVD<Eigen::MatrixXd>::SingularValuesType vSingular = svdA.singularValues();
  // Build a diagonal matrix with the Inverted Singular values
  // The pseudo inverted singular matrix is easy to compute :
  // is formed by replacing every nonzero entry by its reciprocal (inversing).
  Eigen::VectorXd vPseudoInvertedSingular( svdA.matrixV().cols() );

  for (int iRow=0; iRow<vSingular.rows(); iRow++)
  {
    if ( fabs(vSingular(iRow)) <= rcond )
    {
      vPseudoInvertedSingular(iRow)=0.;
    }
    else
      vPseudoInvertedSingular(iRow)=1./vSingular(iRow);
  }
  // A little optimization here
  Eigen::MatrixXd mAdjointU = svdA.matrixU().adjoint().block( 0, 0, vSingular.rows(), svdA.matrixU().adjoint().cols() );
  // Pseudo-Inversion : V * S * U'
  Eigen::MatrixXd a_pinv = (svdA.matrixV() * vPseudoInvertedSingular.asDiagonal()) * mAdjointU;
  if( flip )
  {
    a = a.transpose();
    a_pinv = a_pinv.transpose();
  }
  return a_pinv;
}

IGL_INLINE void igl::n_polyvector_general(const Eigen::MatrixXd &V,
                             const Eigen::MatrixXi &F,
                             const Eigen::VectorXi& b,
                             const Eigen::MatrixXd& bc,
                             const Eigen::VectorXi &I,
                             Eigen::MatrixXd &output)
{
  Eigen::VectorXi isConstrained = Eigen::VectorXi::Constant(F.rows(),0);
  Eigen::MatrixXd cfW = Eigen::MatrixXd::Constant(F.rows(),bc.cols(),0);

  for(unsigned i=0; i<b.size();++i)
  {
    isConstrained(b(i)) = 1;
    cfW.row(b(i)) << bc.row(i);
  }
  int n = I.rows();
  igl::GeneralPolyVectorFieldFinder<Eigen::MatrixXd, Eigen::MatrixXi> pvff(V,F,n);
  pvff.solve(isConstrained, cfW, I, output);
}

//! Function to read the amplitudes and fundamental argument multiplers for tidal corrections
std::pair< Eigen::MatrixXd, Eigen::MatrixXd > readAmplitudesAndFundamentalArgumentMultipliers(
        const std::string amplitudesFile,
        const std::string fundamentalArgumentMultipliersFile,
        const double minimumAmplitude )
{

    // Read amplitudes and fundamental argument multipliers into matrices.
    Eigen::MatrixXd amplitudesRaw = input_output::readMatrixFromFile( amplitudesFile );
    Eigen::MatrixXd fundamentalArgumentMultipliersRaw = input_output::readMatrixFromFile( fundamentalArgumentMultipliersFile );

    // Check whether amplitudes and fundamental argument multipliers matrices have same number of rows
    if( amplitudesRaw.rows( ) != fundamentalArgumentMultipliersRaw.rows( ) )
    {
        throw std::runtime_error( "Amplitude and argument multipler files contain unequal set of entries" );
    }

    // Check whether number of columns in fundamental argument multiplier matrix is equal to 6.
    if( fundamentalArgumentMultipliersRaw.cols( ) != 6 )
    {
        throw std::runtime_error( "Number of columns in fundamental argument multipler matrix not equal to 6" );
    }


    // If acceptance amplitude is larger than 0.0 check RSS of amplitudes for each entry.
    Eigen::MatrixXd fundamentalArgumentMultipliers, amplitudes;
    if( minimumAmplitude > 0.0 )
    {
        // Filter raw data and remove entries with RSS amplitude that is too low
        std::pair< Eigen::MatrixXd, Eigen::MatrixXd > filteredData =
                filterRawDataForAmplitudes( fundamentalArgumentMultipliersRaw, amplitudesRaw, minimumAmplitude );
        amplitudes = filteredData.first;
        fundamentalArgumentMultipliers = filteredData.second;
    }
    // If all amplitudes are to be expected, raw data is same as return data.
    else
    {
        amplitudes = amplitudesRaw;
        fundamentalArgumentMultipliers = fundamentalArgumentMultipliersRaw;
    }

    // return amplitudes and fundamental argument multipliers.
    return std::pair< Eigen::MatrixXd, Eigen::MatrixXd >( amplitudes, fundamentalArgumentMultipliers );
}

void RBMTestCase::inputGradient()
{
    OpenANN::OutputInfo info;
    info.dimensions.push_back(3);
    OpenANN::RBM layer(3, 2, 1, 0.01, 0.0, true);
    LayerAdapter opt(layer, info);

    Eigen::MatrixXd X = Eigen::MatrixXd::Random(2, 3);
    Eigen::MatrixXd Y = Eigen::MatrixXd::Random(2, 2);
    opt.trainingSet(X, Y);
    Eigen::MatrixXd gradient = opt.inputGradient();
    ASSERT_EQUALS(gradient.rows(), 2);
    Eigen::MatrixXd estimatedGradient = OpenANN::FiniteDifferences::
                                        inputGradient(X, Y, opt);
    ASSERT_EQUALS(estimatedGradient.rows(), 2);
    for(int j = 0; j < gradient.rows(); j++)
        for(int i = 0; i < gradient.cols(); i++)
            ASSERT_EQUALS_DELTA(gradient(j, i), estimatedGradient(j, i), 1e-10);
}

//------------------------------------------------------------------------------
void LocalAssembler::assemble_exterior_facet(Eigen::MatrixXd& A,
        UFC& ufc,
        const std::vector<double>& vertex_coordinates,
        const ufc::cell& ufc_cell,
        const Cell& cell,
        const Facet& facet,
        const std::size_t local_facet,
        const MeshFunction<std::size_t>* domains)
{
    // Skip if there are no exterior facet integrals
    if (!ufc.form.has_exterior_facet_integrals())
        return;

    // Extract default exterior facet integral
    ufc::exterior_facet_integral* integral
        = ufc.default_exterior_facet_integral.get();

    // Get integral for sub domain (if any)
    if (domains && !domains->empty())
        integral = ufc.get_exterior_facet_integral((*domains)[facet]);

    // Skip integral if zero
    if (!integral)
        return;

    // Update to current cell
    ufc.update(cell, vertex_coordinates, ufc_cell,
               integral->enabled_coefficients());

    // Tabulate exterior facet tensor
    integral->tabulate_tensor(ufc.A.data(),
                              ufc.w(),
                              vertex_coordinates.data(),
                              local_facet,
                              ufc_cell.orientation);

    // Stuff a_ufc.A into A
    const std::size_t M = A.rows();
    const std::size_t N = A.cols();
    for (std::size_t i = 0; i < M; i++)
        for (std::size_t j = 0; j < N; j++)
            A(i, j) += ufc.A[N*i + j];
}

void BranchList::smoothBranchPositions()
{
	for (int i = 0; i < mBranches.size(); i++)
	{
		Eigen::MatrixXd positions = mBranches[i]->getPositions();
		int numberOfInputPoints = positions.cols();
		int controlPointFactor = 5;
		int numberOfControlPoints = numberOfInputPoints / controlPointFactor;

		vtkCardinalSplinePtr splineX = vtkSmartPointer<vtkCardinalSpline>::New();
		vtkCardinalSplinePtr splineY = vtkSmartPointer<vtkCardinalSpline>::New();
		vtkCardinalSplinePtr splineZ = vtkSmartPointer<vtkCardinalSpline>::New();

		if (numberOfControlPoints >= 2)
		{
			//add control points to spline
			for(int j=0; j<numberOfControlPoints; j++)
			{
				int indexP = (j*numberOfInputPoints)/numberOfControlPoints;

				splineX->AddPoint(indexP,positions(0,indexP));
				splineY->AddPoint(indexP,positions(1,indexP));
				splineZ->AddPoint(indexP,positions(2,indexP));
			}
			//Always add the last point to complete spline
			splineX->AddPoint(numberOfInputPoints-1,positions(0,numberOfInputPoints-1));
			splineY->AddPoint(numberOfInputPoints-1,positions(1,numberOfInputPoints-1));
			splineZ->AddPoint(numberOfInputPoints-1,positions(2,numberOfInputPoints-1));

			//evaluate spline - get smoothed positions
			Eigen::MatrixXd smoothingResult(3 , numberOfInputPoints);
			for(int j=0; j<numberOfInputPoints; j++)
			{
				double splineParameter = j;
				smoothingResult(0,j) = splineX->Evaluate(splineParameter);
				smoothingResult(1,j) = splineY->Evaluate(splineParameter);
				smoothingResult(2,j) = splineZ->Evaluate(splineParameter);
			}
			mBranches[i]->setPositions(smoothingResult);
		}
	}
}

//! Initialize atmosphere table reader.
void TabulatedAtmosphere::initialize( const std::string& atmosphereTableFile )
{
    // Locally store the atmosphere table file name.
    atmosphereTableFile_ = atmosphereTableFile;

    Eigen::MatrixXd containerOfAtmosphereTableFileData
            = input_output::readMatrixFromFile( atmosphereTableFile_, " \t", "%" );

    // Check whether data is present in the file.
    if ( containerOfAtmosphereTableFileData.rows( ) < 1
         || containerOfAtmosphereTableFileData.cols( ) < 1 )
    {
        std::cerr << "The atmosphere table file is empty." << std::endl;
        std::cerr << atmosphereTableFile_ << std::endl;
    }

    // Initialize vectors.
    altitudeData_.resize( containerOfAtmosphereTableFileData.rows( ) );
    densityData_.resize( containerOfAtmosphereTableFileData.rows( ) );
    pressureData_.resize( containerOfAtmosphereTableFileData.rows( ) );
    temperatureData_.resize( containerOfAtmosphereTableFileData.rows( ) );

    // Loop through all the strings stored in the container and store the data
    // in the right Eigen::VectorXd.
    for ( int i = 0; i < containerOfAtmosphereTableFileData.rows( ); i++  )
    {
        altitudeData_[ i ] = containerOfAtmosphereTableFileData( i, 0 );
        densityData_[ i ] = containerOfAtmosphereTableFileData( i, 1 );
        pressureData_[ i ] = containerOfAtmosphereTableFileData( i, 2 );
        temperatureData_[ i ] = containerOfAtmosphereTableFileData( i, 3 );
    }

    using namespace interpolators;

    cubicSplineInterpolationForDensity_
            = boost::make_shared< CubicSplineInterpolatorDouble >( altitudeData_, densityData_ );
    cubicSplineInterpolationForPressure_
            = boost::make_shared< CubicSplineInterpolatorDouble >( altitudeData_, pressureData_ );
    cubicSplineInterpolationForTemperature_
            = boost::make_shared< CubicSplineInterpolatorDouble >(
                altitudeData_, temperatureData_ );
}

bool pre_draw(igl::opengl::glfw::Viewer & viewer)
{
  using namespace Eigen;
  using namespace std;
  if(viewer.core.is_animating)
  {
    // Interpolate pose and identity
    RotationList anim_pose(pose.size());
    for(int e = 0;e<pose.size();e++)
    {
      anim_pose[e] = pose[e].slerp(anim_t,Quaterniond::Identity());
    }
    // Propagate relative rotations via FK to retrieve absolute transformations
    RotationList vQ;
    vector<Vector3d> vT;
    igl::forward_kinematics(C,BE,P,anim_pose,vQ,vT);
    const int dim = C.cols();
    MatrixXd T(BE.rows()*(dim+1),dim);
    for(int e = 0;e<BE.rows();e++)
    {
      Affine3d a = Affine3d::Identity();
      a.translate(vT[e]);
      a.rotate(vQ[e]);
      T.block(e*(dim+1),0,dim+1,dim) =
        a.matrix().transpose().block(0,0,dim+1,dim);
    }
    // Compute deformation via LBS as matrix multiplication
    U = M*T;

    // Also deform skeleton edges
    MatrixXd CT;
    MatrixXi BET;
    igl::deform_skeleton(C,BE,T,CT,BET);

    viewer.data().set_vertices(U);
    viewer.data().set_edges(CT,BET,sea_green);
    viewer.data().compute_normals();
    anim_t += anim_t_dir;
    anim_t_dir *= (anim_t>=1.0 || anim_t<=0.0?-1.0:1.0);
  }
  return false;
}

IGL_INLINE void igl::ViewerData::add_points(const Eigen::MatrixXd& P,  const Eigen::MatrixXd& C)
{
  Eigen::MatrixXd P_temp;

  // If P only has two columns, pad with a column of zeros
  if (P.cols() == 2)
  {
    P_temp = Eigen::MatrixXd::Zero(P.rows(),3);
    P_temp.block(0,0,P.rows(),2) = P;
  }
  else
    P_temp = P;

  int lastid = points.rows();
  points.conservativeResize(points.rows() + P_temp.rows(),6);
  for (unsigned i=0; i<P_temp.rows(); ++i)
    points.row(lastid+i) << P_temp.row(i), i<C.rows() ? C.row(i) : C.row(C.rows()-1);

  dirty |= DIRTY_OVERLAY_POINTS;
}

// Trying to assign from matrix to noOfPoints
void mitk::LibSVMClassifier::ReadXValues(LibSVM::svm_problem * problem, LibSVM::svm_node** xSpace, const Eigen::MatrixXd &X)
{
  auto noOfPoints = static_cast<int>(X.rows());
  auto features = static_cast<int>(X.cols());

  problem->x = static_cast<LibSVM::svm_node **>(malloc(sizeof(LibSVM::svm_node *)  * noOfPoints));
  (*xSpace) = static_cast<LibSVM::svm_node *> (malloc(sizeof(LibSVM::svm_node) * noOfPoints * (features+1)));

  for (int row = 0; row < noOfPoints; ++row)
  {
    for (int col = 0; col < features; ++col)
    {
      (*xSpace)[row*features + col].index = col;
      (*xSpace)[row*features + col].value = X(row,col);
    }
    (*xSpace)[row*features + features].index = -1;

    problem->x[row] = &((*xSpace)[row*features]);
  }
}

bool init_weights(
  const Eigen::MatrixXd & V,
  const Eigen::MatrixXi & F,
  const Eigen::MatrixXd & C,
  const Eigen::MatrixXi & BE,
  Eigen::MatrixXd & W)
{
  using namespace igl;
  using namespace Eigen;
  // Mesh with samples on skeleton
  // New vertices of tet mesh, V prefaces VV
  MatrixXd VV;
  MatrixXi TT,FF,CE;
  VectorXi P;
  if(!mesh_with_skeleton(V,F,C,P,BE,CE,10,VV,TT,FF))
  {
    return false;
  }
  // List of boundary indices (aka fixed value indices into VV)
  VectorXi b;
  // List of boundary conditions of each weight function
  MatrixXd bc;
  if(!boundary_conditions(VV,TT,C,P,BE,CE,b,bc))
  {
    return false;
  }

  // compute BBW
  // Default bbw data and flags
  BBWData bbw_data;
  bbw_data.active_set_params.max_iter = 4;
  // Weights matrix
  if(!bbw(VV,TT,b,bc,bbw_data,W))
  {
    return false;
  }
  // Normalize weights to sum to one
  normalize_row_sums(W,W);
  W.conservativeResize(V.rows(),W.cols());
  return true;
}

double accuracy( const Eigen::MatrixXd& X, const Eigen::VectorXd& y, const Eigen::VectorXd& theta ) {
  const unsigned N = X.rows(), K = X.cols();
  // Eigen:::VectorXd p;
  // p.resize(N, 0); p.setZero();

  double correct = 0.0;
  Eigen::VectorXd xtheta = X*theta;
  Eigen::VectorXd predicted;
  maths::Sigmoid()( predicted, xtheta );

      
  std::cout << predicted << std::endl;
  for ( size_t i = 0; i < N; ++i ) {
    int p = predicted[i]> 0.5 ? 1 : 0;
    // printf("%d, %d, %f\n", )
    if ( p == y(i,0) ) ++correct;
  }

  return correct / N;
  
}