C++ DMatrix类代码示例

 DLLEXPORT void XGBoosterBoostOneIter( void *handle, void *dtrain, 
                             float *grad, float *hess, size_t len, int bst_group ){
     Booster *bst = static_cast<Booster*>(handle);
     DMatrix *dtr = static_cast<DMatrix*>(dtrain);
     bst->CheckInit(); dtr->CheckInit(); 
     bst->BoostOneIter( *dtr, grad, hess, len, bst_group );
 }      

double SQPInternal::quad_form(const std::vector<double>& x, const DMatrix& A){
  // Assert dimensions
  casadi_assert(x.size()==A.size1() && x.size()==A.size2());
  
  // Access the internal data of A
  const std::vector<int> &A_rowind = A.rowind();
  const std::vector<int> &A_col = A.col();
  const std::vector<double> &A_data = A.data();
  
  // Return value
  double ret=0;

  // Loop over the rows of A
  for(int i=0; i<x.size(); ++i){
    // Loop over the nonzeros of A
    for(int el=A_rowind[i]; el<A_rowind[i+1]; ++el){
      // Get column
      int j = A_col[el];
      
      // Add contribution
      ret += x[i]*A_data[el]*x[j];
    }
  }
  
  return ret;
}

 DLLEXPORT void XGBoosterUpdateInteract( void *handle, void *dtrain, const char *action ){
     Booster *bst = static_cast<Booster*>(handle);
     DMatrix *dtr = static_cast<DMatrix*>(dtrain); 
     bst->CheckInit(); dtr->CheckInit(); 
     std::string act( action );
     bst->UpdateInteract( act, *dtr );
 }

//--------------------------------------------------------------------------
double Performance(const DVector &ytest_truth, const DMatrix &Ytest_predict,
                   int NumClasses) {

  long n = ytest_truth.GetN();
  if (n != Ytest_predict.GetM() || Ytest_predict.GetN() != NumClasses ) {
    printf("Performance. Error: Size mismatch. Return NAN");
    return NAN;
  }
  if (NumClasses <= 2) {
    printf("Performance. Error: Not multiclass classification. Return NAN");
    return NAN;
  }
  double perf = 0.0;

  // Compute ytest_predict
  DVector ytest_predict(n);
  double *y = ytest_predict.GetPointer();
  for (long i = 0; i < n; i++) {
    DVector row(NumClasses);
    Ytest_predict.GetRow(i, row);
    long idx = -1;
    row.Max(idx);
    y[i] = (double)idx;
  }

  // Accuracy
  double *y1 = ytest_truth.GetPointer();
  double *y2 = ytest_predict.GetPointer();
  for (long i = 0; i < n; i++) {
    perf += ((int)y1[i])==((int)y2[i]) ? 1.0 : 0.0;
  }
  perf = perf/n * 100.0;

  return perf;
}

returnValue ModelData::setNARXmodel( const uint _delay, const DMatrix& _parms ) {

	if( rhs_name.empty() && NX2 == 0 && NX3 == 0 ) {
		NX2 = _parms.getNumRows();
		delay = _delay;
		uint numParms = 1;
		uint n = delay*getNX();
		if( delay >= 1 ) {
			numParms = numParms + n;
		}
		for( uint i = 1; i < delay; i++ ) {
			numParms = numParms + (n+1)*(uint)pow((double)n,(int)i)/2;
		}
		if( _parms.getNumCols() == numParms ) {
			parms = _parms;
		}
		else {
			return ACADOERROR( RET_UNABLE_TO_EXPORT_CODE );
		}

		export_rhs = BT_TRUE;
	}
	else {
		return ACADOERROR( RET_INVALID_OPTION );
	}
	return SUCCESSFUL_RETURN;
}

returnValue ConstraintElement::get(Function& function_, DMatrix& lb_, DMatrix& ub_)
{
	if ( fcn == NULL )
		return RET_INITIALIZE_FIRST;

	// This is not exactly bullet proof, but for now will serve the purpose
	function_ = fcn[ 0 ];

	int dimFcn = fcn[ 0 ].getDim();

	if ( dimFcn == 0 )
	{
		lb_.init(0, 0);
		ub_.init(0, 0);

		return SUCCESSFUL_RETURN;
	}

	lb_.init(nB, dimFcn);
	ub_.init(nB, dimFcn);

	int i, j;

	for (i = 0; i < nB; ++i)
		for (j = 0; j < dimFcn; ++j)
		{
			lb_(i, j) = lb[ i ][ j ];
			ub_(i, j) = ub[ i ][ j ];
		}

	return SUCCESSFUL_RETURN;
}

void GLEABC::set_BBT(const DMatrix& nBBT)
{
    if (n==0) set_size(nBBT.rows());
    if (A.rows()>0) fr_init=true;
    if (n!=nBBT.rows()) ERROR("Use set_size if you want to change matrix size.");
    if (n!=nBBT.cols()) ERROR("Square matrix expected.");
    BBT=nBBT; fr_c=false;
}

int EigenVectors::TrimSystem(DMatrix &dmat, DArray &darr) const {
    dmat.pop_back();
    darr.pop_back();
    int i, size = dmat.size();
    for(i = 0; i < size; ++i) {
        dmat[i].pop_back();
    }
    return size;
}

returnValue ExportIRK4StageSimplifiedNewton::setEigenvalues( const DMatrix& _eig ) {
    eig = _eig;

    if( _eig.getNumRows() != 2 || _eig.getNumCols() != 2 ) return ACADOERROR( RET_INVALID_ARGUMENTS );
    if( _eig.isZero() ) return ACADOERROR( RET_INVALID_ARGUMENTS );
    // each row represents a complex conjugate pair of eigenvalues

    return SUCCESSFUL_RETURN;
}

void EigenVectors::InitMatrix(DMatrix &dmat, int size) const {
    if(!dmat.empty()) {
        dmat.clear();
    }
    DArray darr(size, 0.);
    dmat.reserve(size);
    for(int i = 0; i < size; ++i) {
        dmat.push_back(darr);
    }
}

//--------------------------------------------------------------------------
void ConvertYtrain(const DVector &ytrain, DMatrix &Ytrain, int NumClasses) {
  Ytrain.Init(ytrain.GetN(), NumClasses);
  // Lingfei: the y label is supposed to start from 0.
  for (int i = 0; i < NumClasses; i++) {
    DVector y = ytrain;
    double *my = y.GetPointer();
    for (long j = 0; j < y.GetN(); j++) {
      my[j] = ((int)my[j])==i ? 1.0 : -1.0;
    }
    Ytrain.SetColumn(i, y);
  }
}

int main(int argc, char *argv[])
{
  /*
  DMatrix W { 5, 5, 
      {   0,   4, inf,   2, inf, 
	inf,   0,   3, inf,   1, 
	  5, inf,   0,   1, inf, 
          7, inf, inf,   0,   5,
	  3, inf, inf, inf,   0 } };
  */

  // Checking that the value of n is given
  if (argc != 2) {
    cerr << "floyd-sequential n" << endl;
    exit(1);
  }
  int n = atoi(argv[1]);
  if (n == 0) {
    cerr << argv[1] << " is a wrong parameter" << endl;
    exit(1);
  }
   
  srand(10);

  DMatrix W = generate_matrix(n);

  DMatrix dist(W.rows(),W.cols(),inf);
  IMatrix next(W.rows(),W.cols(),-1);


  double start = omp_get_wtime();
  floyd_warshall(W, dist, next); 

  
  path_t p = recover_path(dist,next,0,(rand()%(n-1))+1);
  double end = omp_get_wtime();

  if (p.size() == 0) {
    cout << "No path" << endl;
  } else {
    copy(p.begin(), p.end(), ostream_iterator<int>(cout, " "));
    cout << "\n";
    double cost = 0;
    #pragma omp parallel for reduction (+:cost)
    for (uint32_t i = 0; i < p.size()-1; ++i) {
      cost += dist(p[i],p[i+1]);
    }
    cout << "Cost: " << cost << "\n";
  }
  
  cout << "Elapsed time: " << (end-start) << endl;
}

void GLEABC::set_C(const DMatrix& nC)
{
    if (n==0) set_size(nC.rows());
    if (A.rows()>0) fr_init=true;
    if (n!=nC.rows()) ERROR("Use set_size if you want to change matrix size.");
    if (n!=nC.cols()) ERROR("Square matrix expected.");
    C=nC; fr_c=true;
    if (fr_init && fr_c)
    {
        //update BBT, it's so cheap we can do it routinely;
        DMatrix T; mult(A, C, BBT); transpose(BBT, T);  incr(BBT,T);
    }
}

returnValue RungeKuttaExport::initializeButcherTableau( const DMatrix& _AA, const DVector& _bb, const DVector& _cc ) {

	if( _cc.isEmpty() || !_AA.isSquare() || _AA.getNumRows() != _bb.getDim() || _bb.getDim() != _cc.getDim() ) return RET_INVALID_OPTION;

	numStages = _cc.getDim();
	is_symmetric = checkSymmetry( _cc );
//	std::cout << "Symmetry of the chosen method: " << is_symmetric << "\n";
	AA = _AA;
	bb = _bb;
	cc = _cc;

	return SUCCESSFUL_RETURN;
}

  void CSparseCholeskyInternal::prepare() {

    prepared_ = false;

    // Get a reference to the nonzeros of the linear system
    const vector<double>& linsys_nz = input().data();

    // Make sure that all entries of the linear system are valid
    for (int k=0; k<linsys_nz.size(); ++k) {
      casadi_assert_message(!isnan(linsys_nz[k]), "Nonzero " << k << " is not-a-number");
      casadi_assert_message(!isinf(linsys_nz[k]), "Nonzero " << k << " is infinite");
    }

    if (verbose()) {
      userOut() << "CSparseCholeskyInternal::prepare: numeric factorization" << endl;
      userOut() << "linear system to be factorized = " << endl;
      input(0).printSparse();
    }

    if (L_) cs_nfree(L_);
    L_ = cs_chol(&AT_, S_) ;                 // numeric Cholesky factorization
    if (L_==0) {
      DMatrix temp = input();
      temp.makeSparse();
      if (temp.sparsity().issingular()) {
        stringstream ss;
        ss << "CSparseCholeskyInternal::prepare: factorization failed due "
          "to matrix being singular. Matrix contains numerical zeros which are"
          " structurally non-zero. Promoting these zeros to be structural "
          "zeros, the matrix was found to be structurally rank deficient. "
          "sprank: " << sprank(temp.sparsity()) << " <-> " << temp.size2() << endl;
        if (verbose()) {
          ss << "Sparsity of the linear system: " << endl;
          input(LINSOL_A).sparsity().print(ss); // print detailed
        }
        throw CasadiException(ss.str());
      } else {
        stringstream ss;
        ss << "CSparseCholeskyInternal::prepare: factorization failed, "
            "check if Jacobian is singular" << endl;
        if (verbose()) {
          ss << "Sparsity of the linear system: " << endl;
          input(LINSOL_A).sparsity().print(ss); // print detailed
        }
        throw CasadiException(ss.str());
      }
    }
    casadi_assert(L_!=0);

    prepared_ = true;
  }