C++ ColumnVector类代码示例

int LINEAR::recognize( const ColumnVector &v ) const {
  int label;
  if( v.length() < liblinear->prob.n ){
    ERR_PRINT("Dimension of feature vector is too small. %d < %d\n",v.length(),liblinear->prob.n );
    exit(1);
  }
  else if( v.length() > liblinear->prob.n )
    ERR_PRINT("Warning: Dimension of feature vector is too large. %d > %d\n",v.length(),liblinear->prob.n );

  // 特徴ベクトルをセット
  struct feature_node *x = new struct feature_node[ liblinear->prob.n + 2 ];
  int idx = 0;
  for( int i = 0; i < liblinear->prob.n; i++ ){
    x[ idx ].index = i;
    x[ idx ].value = v( i );

    if( is_scaling ){
      // 下記の条件を満たさないときはスキップ（idxを更新しない）	
      if( ( feature_max[ x[ idx ].index ] != feature_min[ x[ idx ].index ] ) && ( x[ idx ].value != 0 ) ){
	x[ idx ].value = scaling( x[ idx ].index, x[ idx ].value );
	if( liblinear->model->bias < 0 ) x[ idx ].index++; // indexが１から始まる
	++idx;
      }
    }
    else{
      // 下記の条件を満たさないときはスキップ（idxを更新しない）	
      if( x[ idx ].value != 0 ){
	if( liblinear->model->bias < 0 ) x[ idx ].index++; // indexが１から始まる
	++idx;
      }
    }
  }

  if(liblinear->model->bias>=0){
    x[ idx ].index = liblinear->prob.n;
    x[ idx ].value = liblinear->model->bias;
    idx++;
  }
  x[ idx ].index = -1;

  // predict
  if( probability ){
    if( liblinear->model->param.solver_type != L2R_LR ){
      ERR_PRINT( "probability output is only supported for logistic regression\n" );
      exit(1);
    }
    label = predict_probability(liblinear->model,x,prob_estimates);
    //     for(j=0;j<nclass;j++)
    //       printf(" %g",prob_estimates[j]);
    //     printf("\n");
  }
  else
    label = predict(liblinear->model,x);

  if( is_y_scaling )
    label = y_scaling( label );

  delete x;  
  return label;
}

void SpinAdapted::xsolve_AxeqB(const Matrix& a, const ColumnVector& b, ColumnVector& x)
{
  FORTINT ar = a.Nrows();
  int bc = 1;
  int info=0;
  FORTINT* ipiv = new FORTINT[ar];
  double* bwork = new double[ar];
  for(int i = 0;i<ar;++i)
    bwork[i] = b.element(i);
  double* workmat = new double[ar*ar];
  for(int i = 0;i<ar;++i)
    for(int j = 0;j<ar;++j)
      workmat[i*ar+j] = a.element(j,i);

  GESV(ar, bc, workmat, ar, ipiv, bwork, ar, info);
  delete[] ipiv;
  delete[] workmat;

  for(int i = 0;i<ar;++i)
    x.element(i) = bwork[i];

  delete[] bwork;

  if(info != 0)
  {
     pout << "Xsolve failed with info error " << info << endl;
     abort();
  }
}

void CActorFromContinuousActionGradientPolicy::receiveError(double critic, CStateCollection *oldState, CAction *Action, CActionData *data)
{
	gradientETraces->updateETraces(Action->getDuration());
	
	CContinuousActionData *contData = NULL;
	if (data)
	{
		contData = dynamic_cast<CContinuousActionData *>(data);
	}
	else
	{
		contData = dynamic_cast<CContinuousActionData *>(Action->getActionData());

	}

	assert(gradientPolicy->getRandomController());
	ColumnVector *noise = gradientPolicy->getRandomController()->getLastNoise();

	if (DebugIsEnabled('a'))
	{
		DebugPrint('a', "ActorCritic Noise: ");
		policyDifference->saveASCII(DebugGetFileHandle('a'));
	}

	for (int i = 0; i < gradientPolicy->getNumOutputs(); i ++)
	{
		gradientFeatureList->clear();
		gradientPolicy->getGradient(oldState, i, gradientFeatureList);
		
		gradientETraces->addGradientETrace(gradientFeatureList, noise->element(i));
	}

	gradientPolicy->updateGradient(gradientETraces->getGradientETraces(), critic * getParameter("ActorLearningRate"));
}

/* Solve Ax = b using the Jacobi Method. */
Matrix jacobi(Matrix& A, Matrix& b)
{
    ColumnVector x0(A.rows()); // our initial guess
    ColumnVector x1(A.rows()); // our next guess
    
    // STEP 1: Choose an initial guess
    fill(x0.begin(),x0.end(),1);
    
    // STEP 2: While convergence is not reached, iterate.
    ColumnVector r = static_cast<ColumnVector>(A*x0 - b);
    while (r.length() > 1)
    {
        for (int i=0;i<A.cols();i++)
        {
            double sum = 0;
            for (int j=0;j<A.cols();j++)
            {
                if (j==i) continue;
                sum = sum + A(i,j) * x0(j,0);
                
            }            
            x1(i,0) = (b(i,0) - sum) / A(i,i);
        }
        x0 = x1;
        r = static_cast<ColumnVector>(A*x0 - b);
    }
    
    shared_ptr<Matrix> final_x(new Matrix(static_cast<Matrix>(x0)));
    return *final_x;
}

CvPoint3D32f Model::getRealCoordinatesFromVoxelMap(int u, int v, int w) {


    //Matrix worldMat(4, 1);
    //worldMat = 0;

    float modelContents[] = { u, v, w, 1 };
    ColumnVector modelMat(4);
    modelMat << modelContents;

    // [x,y,z,1] ?
    // [x,y,z,1] = convMat^{-1} * [u,v,w,1]

    //cvSolve(&A, &b, &x);    // solve (Ax=b) for x
//	cvSolve(convMat, &modelMat, &worldMat);

    ColumnVector worldMat = convMat.i() * modelMat;

    //printCvMat(&worldMat);

    float x = worldMat.element(0);
    float y = worldMat.element(1);
    float z = worldMat.element(2);

    //LOG4CPLUS_DEBUG(myLogger, "(u,v,w)=("<< u <<","<<v<<","<<w<<") -> (x,y,z)=(" << x <<","<<y<<","<<z<<")");

    CvPoint3D32f p = cvPoint3D32f(x, y, z);
    return p;
}

/* solve Ax=b using the Method of Steepest Descent. */
Matrix steepestDescent(Matrix& A, Matrix& b)
{
    // the Method of Steepest Descent *requires* a symmetric matrix.
    if (isSymmetric(A)==false)
    {
        shared_ptr<Matrix> nullMat(new Matrix(0,0));
        return *nullMat;
    }
    
    /* STEP 1: Start with a guess. Our guess is all ones. */
    ColumnVector x(A.cols());
    fill(x.begin(),x.end(),1);
    
    /* This is NOT an infinite loop. There's a break statement inside. */
    while(true)
    {
        /* STEP 2: Calculate the residual r_0 = b - Ax_0 */
        ColumnVector r =  static_cast<ColumnVector> (b - A*x);

        if (r.length() < .01) break;
        
        /* STEP 3: Calculate alpha */
        double alpha = (r.transpose() * r)(0,0) / (r.transpose() * A * r)(0,0);
                
        /* STEP 4: Calculate new X_1 where X_1 = X_0 + alpha*r_0 */
        x = x + alpha * r;
    }
    
    shared_ptr<Matrix> final_x(new Matrix(static_cast<Matrix>(x)));
    
    return *final_x;
}

/* Solve Ax = b using the conjugate gradient method. */
Matrix conjugateGradient(Matrix& A, Matrix& b)
{
    double error_tol = .5;      // error tolerance
    int max_iter = 200;          // max # of iterations
    ColumnVector x(A.rows());   // the solution we will iteratively arrive at
    
    int i = 0;
    ColumnVector r = static_cast<ColumnVector>(b - A*x);
    ColumnVector d = r;
    double sigma_old = 0; // will be used later on, in the loop
    double sigma_new = (r.transpose() * r)(0,0);
    double sigma_0 = sigma_new;
    
    while (i < max_iter && sigma_new > error_tol * error_tol * sigma_0)
    {
        ColumnVector q = A * d;
        double alpha = sigma_new / (d.transpose() * q)(0,0);
        x = x + alpha * d;
        
        if (i % 50 == 0)
        {
            r = static_cast<ColumnVector>(b - A*x);
        }else{
            r = r - alpha * q;
        }
        sigma_old = sigma_new;
        sigma_new = (r.transpose() * r)(0,0);
        double beta = sigma_new / sigma_old;
        d = r + beta * d;
        i++;
    }
    
    shared_ptr<Matrix> final_x(new Matrix(static_cast<Matrix>(x)));    
    return *final_x;
}

// Matrix A's first n columns are orthonormal
// so A.Columns(1,n).t() * A.Columns(1,n) is the identity matrix.
// Fill out the remaining columns of A to make them orthonormal
// so A.t() * A is the identity matrix 
void extend_orthonormal(Matrix& A, int n)
{
   REPORT
   Tracer et("extend_orthonormal");
   int nr = A.nrows(); int nc = A.ncols();
   if (nc > nr) Throw(IncompatibleDimensionsException(A));
   if (n > nc) Throw(IncompatibleDimensionsException(A));
   ColumnVector SSR;
   { Matrix A1 = A.Columns(1,n); SSR = A1.sum_square_rows(); }
   for (int i = n; i < nc; ++i)
   {
      // pick row with smallest SSQ
      int k; SSR.minimum1(k);
      // orthogonalise column with 1 at element k, 0 elsewhere
      // next line is rather inefficient
      ColumnVector X = - A.Columns(1, i) * A.SubMatrix(k, k, 1, i).t();
      X(k) += 1.0;
      // normalise
      X /= sqrt(X.SumSquare());
      // update row sums of squares
      for (k = 1; k <= nr; ++k) SSR(k) += square(X(k));
      // load new column into matrix
      A.Column(i+1) = X;
   }
}

double CLocalRBFRegression::doRegression(ColumnVector *vector, DataSubset *subset)
{
//	cout << "NNs for Input " << vector->t() << endl;
	
	DataSubset::iterator it = subset->begin();
/*	for (int i = 0; it != subset->end(); it ++, i++)
	{
		ColumnVector dist(*vector);
		dist = dist - *(*input)[*it];
		printf("(%d %f) ", *it, dist.norm_Frobenius());
	}
	printf("\n"); */

	ColumnVector *rbfFactors = getRBFFactors(vector, subset);
	it = subset->begin();
	double value = 0;
	for (int i = 0; it != subset->end(); it ++, i++)
	{
		value += rbfFactors->element(i) * (*output)[*it];
	}
	double sum = rbfFactors->sum();
	if (sum > 0)
	{
		value = value / sum;
	}
	//printf("Value: %f %f ", value ,sum);
	//cout << rbfFactors->t();
	return value;
}

ColumnVector find_top_n(ColumnVector W, int n, double cindex){
	ColumnVector out(n);
	int count=0;
	double max=0.0;
	out(0)=cindex;
	for(int i=0;i<W.rows();i++){
		if(W(i)!=0.0) count +=1;
		if(W(i)<0.0) W(i)=(-1)*W(i);
	}
	out(1)=count;
	int k=2;
	while(k<n){
		out(k)=0;
		max=W(0);
		for(int i=1;i<W.rows();i++){
			if(W(i)>max){
				max=W(i);
				out(k)=i;
			}
		}
		W(out(k))=0.0;
		k +=1;
	}
	return out;
}

ColumnVector<Type,N> RowVector<Type,N>::T() const
{
    Trace("RowVector<Type,N>", "T()");

    ColumnVector<Type,N> Result;
    for (unsigned int i = 1; i <= N; ++i) Result.Value(i) = this->Value(i);
    return Result;
}

void test1(Real* y, Real* x1, Real* x2, int nobs, int npred)
{
   cout << "\n\nTest 1 - traditional, bad\n";

   // traditional sum of squares and products method of calculation
   // but not adjusting means; maybe subject to round-off error

   // make matrix of predictor values with 1s into col 1 of matrix
   int npred1 = npred+1;        // number of cols including col of ones.
   Matrix X(nobs,npred1);
   X.Column(1) = 1.0;

   // load x1 and x2 into X
   //    [use << rather than = when loading arrays]
   X.Column(2) << x1;  X.Column(3) << x2;

   // vector of Y values
   ColumnVector Y(nobs); Y << y;

   // form sum of squares and product matrix
   //    [use << rather than = for copying Matrix into SymmetricMatrix]
   SymmetricMatrix SSQ; SSQ << X.t() * X;

   // calculate estimate
   //    [bracket last two terms to force this multiplication first]
   //    [ .i() means inverse, but inverse is not explicity calculated]
   ColumnVector A = SSQ.i() * (X.t() * Y);

   // Get variances of estimates from diagonal elements of inverse of SSQ
   // get inverse of SSQ - we need it for finding D
   DiagonalMatrix D; D << SSQ.i();
   ColumnVector V = D.AsColumn();

   // Calculate fitted values and residuals
   ColumnVector Fitted = X * A;
   ColumnVector Residual = Y - Fitted;
   Real ResVar = Residual.SumSquare() / (nobs-npred1);

   // Get diagonals of Hat matrix (an expensive way of doing this)
   DiagonalMatrix Hat;  Hat << X * (X.t() * X).i() * X.t();

   // print out answers
   cout << "\nEstimates and their standard errors\n\n";

   // make vector of standard errors
   ColumnVector SE(npred1);
   for (int i=1; i<=npred1; i++) SE(i) = sqrt(V(i)*ResVar);
   // use concatenation function to form matrix and use matrix print
   // to get two columns
   cout << setw(11) << setprecision(5) << (A | SE) << endl;

   cout << "\nObservations, fitted value, residual value, hat value\n";

   // use concatenation again; select only columns 2 to 3 of X
   cout << setw(9) << setprecision(3) <<
     (X.Columns(2,3) | Y | Fitted | Residual | Hat.AsColumn());
   cout << "\n\n";
}

void NonLinearLeastSquares::Fit(const ColumnVector& Data,
   ColumnVector& Parameters)
{
   Tracer tr("NonLinearLeastSquares::Fit");
   n_param = Parameters.Nrows(); n_obs = Data.Nrows();
   DataPointer = &Data;
   FindMaximum2::Fit(Parameters, Lim);
   cout << "\nConverged\n";
}

/** 
 * Compute a distance metric between two columns of a
 * matrix. <b>Note that the indexes are *1* based (not 0) as that is
 * Newmat's convention</b>. Note that dist(M,i,j) must equal dist(M,j,i);
 * 
 * @param M - Matrix whose columns represent individual items to be clustered.
 * @param col1Ix - Column index to be compared (1 based).
 * @param col2Ix - Column index to be compared (1 based).
 * 
 * @return - "Distance" or "dissimilarity" metric between two columns of matrix.
 */
double GuassianRadial::dist(const Matrix &M, int col1Ix, int col2Ix) const {
  double dist = 0;
  if(col1Ix == col2Ix) 
    return 0;
  ColumnVector V = M.Column(col1Ix) - M.Column(col2Ix);
  dist = V.SumSquare() / (2 * m_Sigma * m_Sigma);
  dist = exp(-1 * dist);
  return dist;
}

const ColumnVector* RowVector::transpose(const RowVector *matA) {
    ColumnVector* t = new ColumnVector(matA->cols());
    
    for (int i = 0; i < matA->cols(); i++) {
        t->set(i, matA->get(i));
    }
    
    return t;
}