C++ ColumnVector::Nrows方法代码示例

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" << endl;
}

std::vector<float> BinghamThread::fit_bingham( const ColumnVector& sh_data,
                                                  const Matrix& tess,
                                                  const std::vector<QSet<int> >& adj,
                                                  const Matrix& base,
                                                  const int neighborhood,
                                                  const int num_max )
{
    unsigned int mod = 9;
    // reserve memory:
    std::vector<float> result( 27, 0 );


    // if no CSD no fit necessary.
    if ( sh_data( 1 ) == 0 )
    {
        return result;
    }

    // get maxima:
    ColumnVector radius = base * sh_data;

    std::vector<float> qfRadius( radius.Nrows() );
    for ( unsigned int i = 0; i < qfRadius.size(); ++i )
    {
        qfRadius[i] = radius( i + 1 );
    }

    std::vector<int> qiRadius( radius.Nrows() );
    for ( unsigned int i = 0; i < qiRadius.size(); ++i )
    {
        qiRadius[i] = i;
    }

    std::vector<int> maxima;
    for ( unsigned int i = 0; i < qfRadius.size(); ++i )
    {
        QSet<int> n = adj[i];
        float r = qfRadius[i];
        if ( r > 0 )
        {
            bool isMax = true;
            foreach (const int &value, n)
            {
                if ( r < qfRadius[value] )
                {
                    isMax = false;
                }
            }
            if ( isMax )
            {
                maxima.push_back( i );
            }
        }
    }

void SpectClust::multByMatrix(ColumnVector &N, ColumnVector &O, const Matrix &M) {
  int nRow = M.Nrows(), nCol = M.Ncols();
  if(!(N.Nrows() == O.Nrows() && M.Nrows() == M.Ncols() && M.Nrows() == N.Nrows())) {
    Err::errAbort("wrong dimensions: " + ToStr(O.Nrows()) + " " + ToStr(N.Nrows()) + " " + ToStr(M.Nrows()));
  }
  for(int rowIx = 0; rowIx < nRow; rowIx++) {
    N[rowIx] = 0;
    for(int colIx = 0; colIx < nCol; colIx++) {
      N[rowIx] += O[colIx] * M[rowIx][colIx];
    }
  }
}

void FFT(const ColumnVector& U, const ColumnVector& V,
   ColumnVector& X, ColumnVector& Y)
{
   // from Carl de Boor (1980), Siam J Sci Stat Comput, 1 173-8
   // but first try Sande and Gentleman
   Tracer trace("FFT");
   REPORT
   const int n = U.Nrows();                     // length of arrays
   if (n != V.Nrows() || n == 0)
      Throw(ProgramException("Vector lengths unequal or zero", U, V));
   if (n == 1) { REPORT X = U; Y = V; return; }

   // see if we can use the newfft routine
   if (!FFT_Controller::OnlyOldFFT && FFT_Controller::CanFactor(n))
   {
      REPORT
      X = U; Y = V;
      if ( FFT_Controller::ar_1d_ft(n,X.Store(),Y.Store()) ) return;
   }

   ColumnVector B = V;
   ColumnVector A = U;
   X.ReSize(n); Y.ReSize(n);
   const int nextmx = 8;
#ifndef ATandT
   int prime[8] = { 2,3,5,7,11,13,17,19 };
#else
   int prime[8];
   prime[0]=2; prime[1]=3; prime[2]=5; prime[3]=7;
   prime[4]=11; prime[5]=13; prime[6]=17; prime[7]=19;
#endif
   int after = 1; int before = n; int next = 0; bool inzee = true;
   int now = 0; int b1;             // initialised to keep gnu happy

   do
   {
      for (;;)
      {
     if (next < nextmx) { REPORT now = prime[next]; }
     b1 = before / now;  if (b1 * now == before) { REPORT break; }
     next++; now += 2;
      }
      before = b1;

      if (inzee) { REPORT fftstep(A, B, X, Y, after, now, before); }
      else { REPORT fftstep(X, Y, A, B, after, now, before); }

      inzee = !inzee; after *= now;
   }

void RealFFT(const ColumnVector& U, ColumnVector& X, ColumnVector& Y)
{
   // Fourier transform of a real series
   Tracer trace("RealFFT");
   REPORT
   const int n = U.Nrows();                     // length of arrays
   const int n2 = n / 2;
   if (n != 2 * n2)
      Throw(ProgramException("Vector length not multiple of 2", U));
   ColumnVector A(n2), B(n2);
   Real* a = A.Store(); Real* b = B.Store(); Real* u = U.Store(); int i = n2;
   while (i--) { *a++ = *u++; *b++ = *u++; }
   FFT(A,B,A,B);
   int n21 = n2 + 1;
   X.ReSize(n21); Y.ReSize(n21);
   i = n2 - 1;
   a = A.Store(); b = B.Store();              // first els of A and B
   Real* an = a + i; Real* bn = b + i;        // last els of A and B
   Real* x = X.Store(); Real* y = Y.Store();  // first els of X and Y
   Real* xn = x + n2; Real* yn = y + n2;      // last els of X and Y

   *x++ = *a + *b; *y++ = 0.0;                // first complex element
   *xn-- = *a++ - *b++; *yn-- = 0.0;          // last complex element

   int j = -1; i = n2/2;
   while (i--)
   {
      Real c,s; cossin(j--,n,c,s);
      Real am = *a - *an; Real ap = *a++ + *an--;
      Real bm = *b - *bn; Real bp = *b++ + *bn--;
      Real samcbp = s * am + c * bp; Real sbpcam = s * bp - c * am;
      *x++  =  0.5 * ( ap + samcbp); *y++  =  0.5 * ( bm + sbpcam);
      *xn-- =  0.5 * ( ap - samcbp); *yn-- =  0.5 * (-bm + sbpcam);
   }
}

void basisfield::SetCoef(const ColumnVector& pcoef) 
{
  if (pcoef.Nrows() != int(CoefSz())) {throw BasisfieldException("basisfield::SetCoef::Mismatch between input vector and # of coefficients");}
  if (!coef) {coef = boost::shared_ptr<NEWMAT::ColumnVector>(new NEWMAT::ColumnVector(pcoef));}
  else {*coef = pcoef;}
  futd.assign(4,false);
}

void FFTI(const ColumnVector& U, const ColumnVector& V,
   ColumnVector& X, ColumnVector& Y)
{
   // Inverse transform
   Tracer trace("FFTI");
   REPORT
   FFT(U,-V,X,Y);
   const Real n = X.Nrows(); X /= n; Y /= (-n);
}

/** Utility printing function. */
void printColumnVector(ColumnVector &v, std::ostream *out, const std::string& delim) {
  int nRow = v.Nrows();
  int i = 0;
  if(out == NULL) 
    out = &cout;
  for(i = 0; i < nRow - 1; i++)
    (*out) << v.element(i) << delim;
  (*out) << v.element(i);
}

double unsupervised::gaussian(int k, const ColumnVector& ob){
	ColumnVector tmp(Dim);
	tmp = ob-mu[k];
	double x = (tmp.t() * sigma[k].i() * tmp).AsScalar();
	if( x<0.0 || std::isnan(x) ) x = 0.0; //for avoiding the failure from calculation error of newmat library.
    
	tmp.CleanUp();
	return (exp(x*(-0.5)) / (pow(2.0*M_PI,ob.Nrows()*0.5)*pow(sigma[k].Determinant(),0.5)));
}

void basisfield::Set(const ColumnVector& pfield) 
{
  if (pfield.Nrows() != int(FieldSz())) {throw BasisfieldException("basisfield::Set::Mismatch between input vector and size of field");}

  volume<float>  vol_pfield(FieldSz_x(),FieldSz_y(),FieldSz_z());
  vol_pfield.setdims(Vxs_x(),Vxs_y(),Vxs_z());
  vol_pfield.insert_vec(pfield);
  
  Set(vol_pfield);
}

//Delete a BV.  Very messy
void SOGP::delete_bv(int loc){
  //First swap loc to the last spot
  RowVector alphastar = alpha.Row(loc);
  alpha.Row(loc)=alpha.Row(alpha.Nrows());
  //Now C
  double cstar = C(loc,loc);
  ColumnVector Cstar = C.Column(loc);
  Cstar(loc)=Cstar(Cstar.Nrows());
  Cstar=Cstar.Rows(1,Cstar.Nrows()-1);
  ColumnVector Crep=C.Column(C.Ncols());
  Crep(loc)=Crep(Crep.Nrows());
  C.Row(loc)=Crep.t();;
  C.Column(loc)=Crep;
  //and Q
  double qstar = Q(loc,loc);
  ColumnVector Qstar = Q.Column(loc);
  Qstar(loc)=Qstar(Qstar.Nrows());
  Qstar=Qstar.Rows(1,Qstar.Nrows()-1);
  ColumnVector Qrep=Q.Column(Q.Ncols());
  Qrep(loc)=Qrep(Qrep.Nrows());
  Q.Row(loc)=Qrep.t();
  Q.Column(loc)=Qrep;

  //Ok, now do the actual removal  Appendix G section g
  alpha= alpha.Rows(1,alpha.Nrows()-1);
  ColumnVector qc = (Qstar+Cstar)/(qstar+cstar);
  for(int i=1;i<=alpha.Ncols();i++)
    alpha.Column(i)-=alphastar(i)*qc;
  C = C.SymSubMatrix(1,C.Ncols()-1) + (Qstar*Qstar.t())/qstar - ((Qstar+Cstar)*(Qstar+Cstar).t())/(qstar+cstar);
  Q = Q.SymSubMatrix(1,Q.Ncols()-1) - (Qstar*Qstar.t())/qstar;
  
  //And the BV
  BV.Column(loc)=BV.Column(BV.Ncols());
  BV=BV.Columns(1,BV.Ncols()-1);
  
  current_size--;
}

ReturnMatrix trans(const ColumnVector & a)
//!  @brief Translation.
{
   Matrix translation(4,4);

   translation << fourbyfourident; // identity matrix 

   if (a.Nrows() == 3) 
     {
       translation(1,4) = a(1);
       translation(2,4) = a(2);
       translation(3,4) = a(3);
     }
   else
       cerr << "trans: wrong size in input vector." << endl;

   translation.Release(); return translation;
}

static void SlowFT(const ColumnVector& a, const ColumnVector&b,
   ColumnVector& x, ColumnVector& y)
{
   int n = a.Nrows();
   x.ReSize(n); y.ReSize(n);
   Real f = 6.2831853071795864769/n;
   for (int j=1; j<=n; j++)
   {
      Real sumx = 0.0; Real sumy = 0.0;
      for (int k=1; k<=n; k++)
      {
	 Real theta = - (j-1) * (k-1) * f;
	 Real c = cos(theta); Real s = sin(theta);
	 sumx += c * a(k) - s * b(k); sumy += s * a(k) + c * b(k);
      }
      x(j) = sumx; y(j) = sumy;
   }
}

static void SlowDTT_II(const ColumnVector& a, ColumnVector& c, ColumnVector& s)
{
   int n = a.Nrows(); c.ReSize(n); s.ReSize(n);
   Real f = 6.2831853071795864769 / (4*n);
   int k;

   for (k=1; k<=n; k++)
   {
      Real sum = 0.0;
      const int k1 = k-1;              // otherwise Visual C++ 5 fails
      for (int j=1; j<=n; j++) sum += cos(k1 * (2*j-1) * f) * a(j);
      c(k) = sum;
   }

   for (k=1; k<=n; k++)
   {
      Real sum = 0.0;
      for (int j=1; j<=n; j++) sum += sin(k * (2*j-1) * f) * a(j);
      s(k) = sum;
   }
}

static void SlowDTT(const ColumnVector& a, ColumnVector& c, ColumnVector& s)
{
   int n1 = a.Nrows(); int n = n1 - 1;
   c.ReSize(n1); s.ReSize(n1);
   Real f = 6.2831853071795864769 / (2*n);
   int k;

   int sign = 1;
   for (k=1; k<=n1; k++)
   {
      Real sum = 0.0;
      for (int j=2; j<=n; j++) sum += cos((j-1) * (k-1) * f) * a(j);
      c(k) = sum + (a(1) + sign * a(n1)) / 2.0;
      sign = -sign;
   }

   for (k=2; k<=n; k++)
   {
      Real sum = 0.0;
      for (int j=2; j<=n; j++) sum += sin((j-1) * (k-1) * f) * a(j);
      s(k) = sum;
   }
   s(1) = s(n1) = 0;
}