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

double getvlog(const Matrix &W, const Matrix &T, const ColumnVector &bf,
	       double cterm, int df, double tol)
{
  Matrix Vk, Vtr; ColumnVector gbeta; double lgbeta, ldetVtr, xbfVtri;
  Vk=T.i()*W; gbeta=Vk.i()*bf; QRD Vtrqr(Vk.t(),tol); Vtr=Vtrqr.R();
  lgbeta=log(fabs(gbeta(1))); //Rprintf("log(abs(gbeta[1])): %f\n", lgbeta);
  ldetVtr=log(fabs(Vtr.Determinant())); //Rprintf("ldetVtr: %f\n",ldetVtr);
  xbfVtri=(bf.t()*Vtr.i()*(bf.t()*Vtr.i()).t()).AsScalar(); //Rprintf("xbfVtri: %f\n",xbfVtri);
  return cterm-ldetVtr+df*lgbeta-0.5*df*xbfVtri;
}

double SpinAdapted::dotproduct(const ColumnVector& a, const ColumnVector& b)
{
  assert(a.Nrows() == b.Nrows());
#ifdef BLAS
  return DDOT(a.Storage(), a.Store(), 1, b.Store(), 1);
#else
  return a.t() * b;
#endif
}

bool SpectClust::MaxEigen(const Matrix &M, double &maxValue, ColumnVector &MaxVec, int maxIterations) {
  double maxDelta = 1e-6;
  bool converged = false;
  int i = 0;
  int nRows = M.Ncols();
  if(M.Ncols() != M.Nrows()) 
    Err::errAbort("MaxEigen() - Can't get eigen values of non square matrices.");
  if(M.Ncols() <= 0) 
    Err::errAbort("MaxEigen() - Must have positive number of rows and columns.");
  ColumnVector V(M.Ncols());
  V = 1.0 / M.Nrows(); // any vector really...
  V = V / Norm1(V);
  MaxVec.ReSize(M.Ncols());
  for(i = 0; i < maxIterations; ++i) {
    //    MaxVec = M * V;
    multByMatrix(MaxVec, V, M);
    double delta = 0;
    double norm = sqrt(SumSquare(MaxVec));
    for(int vIx = 0; vIx < nRows; vIx++) {
      MaxVec.element(vIx) = MaxVec.element(vIx) / norm; // scale so we don't get too big.
    }
    for(int rowIx = 0; rowIx < nRows; rowIx++) {
      delta += fabs((double)MaxVec.element(rowIx) - V.element(rowIx));
    }
    if(delta < maxDelta)  
      break; // we've already converged to eigen vector.
    V = MaxVec;
  }
  if(i < maxIterations) {
    converged = true;
  }
  // calculate approximate max eigen value using Rayleigh quotient (x'*M*x/x'*x).
  Matrix num = (MaxVec.t() * M * MaxVec);
  Matrix denom =  (MaxVec.t() * MaxVec);
  maxValue = num.element(0,0) / denom.element(0,0);
  return converged;
}

void CFeasibilityMap::initilize_D_and_S(CData &Data) {
	
  ColumnVector order_to_test = get_order_to_test(Data.n_tau, Data.n_var);
  ColumnVector list_feasible_type2 = get_feasible_tau(Data);
	
	for (int i_faulty = 1; i_faulty <= Data.n_faulty; i_faulty++) {  
		
    bool is_pass = false;
    int i_original = Data.Faulty2Original[i_faulty-1];
    ColumnVector x_tilde_i = (Data.D_Observed.row(i_original)).t(); 
    for (int i_order = 1; i_order <= order_to_test.nrows() && !is_pass; i_order++) {
      int i_tau = order_to_test(i_order);
      ColumnVector s_i = tau_to_s_fn(i_tau,Data.n_var);
    
      bool skip_for_type2 = false; 
      if (Data.is_case(i_original,2) && list_feasible_type2(i_tau) == 0) { skip_for_type2 = true;}
      if (!skip_for_type2){
        ColumnVector x_mean;
        int is_feasible = feasible_test_fn(Data, x_tilde_i, s_i, i_original, true, Data.epsilon, x_mean);
        
        if (is_feasible > 0) {
          //copy solution to a temp vector
          ColumnVector temp = x_tilde_i;
          for (int index = 1, count =0; index <= Data.n_var; index++){
            if (s_i(index) == 1) {
              count++;
              temp(index) = x_mean(count);
            }
          }
          Data.Debug = Debug;
          if (Data.PassEdits(temp)) { //then check if it satifies edits          
            is_pass = true;
            Data.initial_S_Mat.row(i_faulty) = s_i.t();   
            Data.D_initial.row(i_original) = temp.t();
          }
        }//is_feasible
        Debug = false;
      }
    } 
  
	}  
	
}

bool SpectClust::findNLargestEvals(const Matrix &M, int numLamda, std::vector<Numeric> &eVals, Matrix &EVec, int maxIterations) {
  bool converged = true;
  EVec.ReSize(M.Ncols(), numLamda);
  eVals.clear();
  eVals.reserve(numLamda);
  Matrix W = M;
    
  for(int i = 1; i <= numLamda; i++) {
    ColumnVector maxVec;
    double maxVal;
    /* Get the maximum eigen vector. */
    converged = MaxEigen(W, maxVal, maxVec, maxIterations) && converged;
    EVec.Column(i) << maxVec;
    eVals.push_back(maxVal);
     
    /* Now subtract of the largest eigen value to get the next
       largest in next iteration. */
    Matrix ToSub = maxVal * (maxVec * maxVec.t());
    W = W - ToSub;
  }
  return converged;
}

//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--;
}

//Add this input and output to the GP
void SOGP::add(const ColumnVector& in,const ColumnVector& out){
  double kstar =m_params.m_kernel->kstar(in);
  
  if(current_size==0){//First point is easy
    C.ReSize(1,1);
    Q.ReSize(1,1);
    //Equations 2.46 with q, r, and s collapsed
    alpha=out.t()/(kstar+m_params.s20);
    C(1,1)=-1/(kstar+m_params.s20);
    Q(1,1)=1/kstar;
    
    current_size=1;
    BV = in;
  }
  else{   //We already have data
    //perform the kernel
    ColumnVector k=m_params.m_kernel->kernelM(in,BV);

    RowVector m=k.t()*alpha;
    double s2 = kstar+(k.t()*C*k).AsScalar();
    
    if(s2<1e-12){//For numerical stability..from Csato's Matlab code?
      //printf("SOGP::Small s2 %lf\n",s2);
      s2=1e-12;
    }
    
    //Update scalars
    //page 33 - Assumes Gaussian noise
    double r = -1/(m_params.s20+s2);
    //printf("out %d m %d r %f\n",out.Nrows(),m.Ncols(),r);
    RowVector q = -r*(out.t()-m);

    //projection onto current BV
    ColumnVector ehat = Q*k;//Appendix G, section c
    //residual length
    double gamma = kstar-(k.t()*ehat).AsScalar();//Ibid
    if(gamma<1e-12){//Numerical instability?
      //printf("SOGP::Gamma (%lf) < 0\n",gamma);
      gamma=0;
    }

    if(gamma<1e-6 && m_params.capacity!= -1){//Nearly singular, do a sparse update (e_tol)
      //printf("SOGP::Sparse! %lf \n",gamma);
      double eta = 1/(1+gamma*r);//Ibid
      ColumnVector shat = C*k+ehat;//Appendix G section e
      alpha=alpha+shat*(q*eta);//Appendix G section f
      C=C+r*eta*shat*shat.t();//Ibid
    }
    else{//Full update
      //printf("SOGP::Full!\n");
      //Expansions are messy
      RowVector expander(1);

      //s is of length N+1
      expander(1)=1;
      ColumnVector s = C*k & expander;//Apendix G section e

      //Add a row to alpha
      expander.ReSize(alpha.Ncols());
      for(int i=0;i<alpha.Ncols();i++)
	expander(i+1)=0;
      alpha = alpha & expander;
      //Update alpha
      alpha=alpha+s*q;//Equations 2.46

      //Add Row to C
      expander.ReSize(C.Ncols());
      for(int i=0;i<C.Ncols();i++)
	expander(i+1)=0;
      C = C&expander;
      //Add a Column to C
      expander.ReSize(C.Nrows());
      for(int i=0;i<C.Nrows();i++)
	expander(i+1)=0;
      C = C | expander.t();
      //Update C
      C = C + r*s*s.t();//Ibid

      //Save the data, N++
      BV = BV | in;
      current_size++;
  
      //Add row to Gram Matrix
      expander.ReSize(Q.Ncols());
      for(int i=0;i<Q.Ncols();i++)
	expander(i+1)=0;
      Q = Q&expander;
      //Add column
      expander.ReSize(Q.Nrows());
      for(int i=0;i<Q.Nrows();i++)
	expander(i+1)=0;
      Q = Q | expander.t();
      //Add one more to ehat
      expander.ReSize(1);
      expander(1)=-1;
      ehat = ehat & expander;
      //Update gram matrix
      Q = Q + (1/gamma)*ehat*ehat.t();//Equation 3.5
    }
//.........这里部分代码省略.........

void trymat6()
{
   Tracer et("Sixth test of Matrix package");
   Tracer::PrintTrace();

   int i,j;


   DiagonalMatrix D(6);
   UpperTriangularMatrix U(6);
   for (i=1;i<=6;i++) { for (j=i;j<=6;j++) U(i,j)=i*i*i-50; D(i,i)=i*i+i-10; }
   LowerTriangularMatrix L=(U*3.0).t();
   SymmetricMatrix S(6);
   for (i=1;i<=6;i++) for (j=i;j<=6;j++) S(i,j)=i*i+2.0+j;
   Matrix MD=D; Matrix ML=L; Matrix MU=U; Matrix MS=S;
   Matrix M(6,6);
   for (i=1;i<=6;i++) for (j=1;j<=6;j++) M(i,j)=i*j+i*i-10.0;  
   {
      Tracer et1("Stage 1");
      Print(Matrix(MS+(-MS)));
      Print(Matrix((S+M)-(MS+M)));
      Print(Matrix((M+U)-(M+MU)));
      Print(Matrix((M+L)-(M+ML)));
   }
   {
      Tracer et1("Stage 2");
      Print(Matrix((M+D)-(M+MD)));
      Print(Matrix((U+D)-(MU+MD)));
      Print(Matrix((D+L)-(ML+MD)));
      Print(Matrix((-U+D)+MU-MD));
      Print(Matrix((-L+D)+ML-MD));
   }
   {
      Tracer et1("Stage 3 - concatenate");
      RowVector A(5);
      A << 1 << 2 << 3 << 4 << 5;
      RowVector B(5);
      B << 3 << 1 << 4 << 1 << 5;
      Matrix C(3,5);
      C <<  2 <<  3 <<  5 <<  7 << 11
        << 13 << 17 << 19 << 23 << 29
        << 31 << 37 << 41 << 43 << 47;
      Matrix X1 = A & B & C;
      Matrix X2 = (A.t() | B.t() | C.t()).t();
      Matrix X3(5,5);
      X3.Row(1)=A; X3.Row(2)=B; X3.Rows(3,5)=C;
      Print(Matrix(X1-X2));
      Print(Matrix(X1-X3));
      LowerTriangularMatrix LT1; LT1 << (A & B & C);
      UpperTriangularMatrix UT1; UT1 << (A.t() | B.t() | C.t());
      Print(LowerTriangularMatrix(LT1-UT1.t()));
      DiagonalMatrix D1; D1 << (A.t() | B.t() | C.t());
      ColumnVector At = A.t();
      ColumnVector Bt = B.t();
      Matrix Ct = C.t();
      LowerTriangularMatrix LT2; LT2 << (At | Bt | Ct);
      UpperTriangularMatrix UT2; UT2 << (At.t() & Bt.t() & Ct.t());
      Matrix ABt = At | Bt;
      DiagonalMatrix D2; D2 << (ABt | Ct);
      Print(LowerTriangularMatrix(LT2-UT2.t()));
      Print(DiagonalMatrix(D1-D2));
      Print(Matrix(LT1+UT2-D2-X1));
      Matrix M1 = LT1 | UT2; Matrix M2 = UT1 & LT2;
      Print(Matrix(M1-M2.t()));
      M1 = UT2 | LT1; M2 = LT2 & UT1;
      Print(Matrix(M1-M2.t()));
      M1 = (LT1 | UT2) & (UT2 | LT1);
      M2 = (UT1 & LT2) | (LT2 & UT1);
      Print(Matrix(M1-M2.t()));
      SymmetricMatrix SM1; SM1 << (M1 + M1.t());
      SymmetricMatrix SM2; SM2 << ((SM1 | M1) & (M1.t() | SM1));
      Matrix M3(20,20);
      M3.SubMatrix(1,10,1,10) = SM1;
      M3.SubMatrix(1,10,11,20) = M1;
      M3.SubMatrix(11,20,1,10) = M2;
      M3.SubMatrix(11,20,11,20) = SM1;
      Print(Matrix(M3-SM2));

      SymmetricMatrix SM(15); SM = 0; SM.SymSubMatrix(1,10) = SM1;
      M3.ReSize(15,15); M3 = 0; M3.SubMatrix(1,10,1,10) = SM1;
      M3 -= SM; Print(M3);
      SM = 0; SM.SymSubMatrix(6,15) = SM1;
      M3.ReSize(15,15); M3 = 0; M3.SubMatrix(6,15,6,15) = SM1;
      M3 = M3.t() - SM; Print(M3);
   }
   {
      Tracer et1("Stage 4 - sort");
      TestSort(1); TestSort(2); TestSort(3); TestSort(4);
      TestSort(15); TestSort(16); TestSort(17); TestSort(18);
      TestSort(99); TestSort(100); TestSort(101);
   }


//   cout << "\nEnd of sixth test\n";
}

void trymat2()
{
//   cout << "\nSecond test of Matrix package\n\n";
   Tracer et("Second test of Matrix package");
   Tracer::PrintTrace();

   int i,j;

   Matrix M(3,5);
   for (i=1; i<=3; i++) for (j=1; j<=5; j++) M(i,j) = 100*i + j;
   Matrix X(8,10);
   for (i=1; i<=8; i++) for (j=1; j<=10; j++) X(i,j) = 1000*i + 10*j;
   Matrix Y = X; Matrix Z = X;
   { X.SubMatrix(2,4,3,7) << M; }
   for (i=1; i<=3; i++) for (j=1; j<=5; j++) Y(i+1,j+2) = 100*i + j;
   Print(Matrix(X-Y));


   Real a[15]; Real* r = a;
   for (i=1; i<=3; i++) for (j=1; j<=5; j++) *r++ = 100*i + j;
   { Z.SubMatrix(2,4,3,7) << a; }
   Print(Matrix(Z-Y));

   { M=33; X.SubMatrix(2,4,3,7) << M; }
   { Z.SubMatrix(2,4,3,7) = 33; }
   Print(Matrix(Z-X));

   for (i=1; i<=8; i++) for (j=1; j<=10; j++) X(i,j) = 1000*i + 10*j;
   Y = X;
   UpperTriangularMatrix U(5);
   for (i=1; i<=5; i++) for (j=i; j<=5; j++) U(i,j) = 100*i + j;
   { X.SubMatrix(3,7,5,9) << U; }
   for (i=1; i<=5; i++) for (j=i; j<=5; j++) Y(i+2,j+4) = 100*i + j;
   for (i=1; i<=5; i++) for (j=1; j<i; j++) Y(i+2,j+4) = 0.0;
   Print(Matrix(X-Y));
   for (i=1; i<=8; i++) for (j=1; j<=10; j++) X(i,j) = 1000*i + 10*j;
   Y = X;
   for (i=1; i<=5; i++) for (j=i; j<=5; j++) U(i,j) = 100*i + j;
   { X.SubMatrix(3,7,5,9).Inject(U); }
   for (i=1; i<=5; i++) for (j=i; j<=5; j++) Y(i+2,j+4) = 100*i + j;
   Print(Matrix(X-Y));


   // test growing and shrinking a vector
   {
      ColumnVector V(100);
      for (i=1;i<=100;i++) V(i) = i*i+i;
      V = V.Rows(1,50);               // to get first 50 vlaues.

      {
         V.Release(); ColumnVector VX=V;
         V.ReSize(100); V = 0.0; V.Rows(1,50)=VX;
      }                               // V now length 100

      M=V; M=100;                     // to make sure V will hold its values
      for (i=1;i<=50;i++) V(i) -= i*i+i;
      Print(V);


	   // test redimensioning vectors with two dimensions given
      ColumnVector CV1(10); CV1 = 10;
      ColumnVector CV2(5); CV2.ReSize(10,1); CV2 = 10;
      V = CV1-CV2; Print(V);

      RowVector RV1(20); RV1 = 100;
      RowVector RV2; RV2.ReSize(1,20); RV2 = 100;
      V = (RV1-RV2).t(); Print(V);

      X.ReSize(4,7);
      for (i=1; i<=4; i++) for (j=1; j<=7; j++) X(i,j) = 1000*i + 10*j;
      Y = 10.5 * X;
      Z = 7.25 - Y;
      M = Z + X * 10.5 - 7.25;
      Print(M);
      Y = 2.5 * X;
      Z = 9.25 + Y;
      M = Z - X * 2.5 - 9.25;
      Print(M);
      U.ReSize(8);
      for (i=1; i<=8; i++) for (j=i; j<=8; j++) U(i,j) = 100*i + j;
      Y = 100 - U;
      M = Y + U - 100;
      Print(M);
   }

   {
      SymmetricMatrix S,T;

      S << (U + U.t());
      T = 100 - S; M = T + S - 100; Print(M);
      T = 100 - 2 * S; M = T + S * 2 - 100; Print(M);
      X = 100 - 2 * S; M = X + S * 2 - 100; Print(M);
      T = S; T = 100 - T; M = T + S - 100; Print(M);
   }

   // test new
   {
      ColumnVector CV1; RowVector RV1;
      Matrix* MX; MX = new Matrix; if (!MX) Throw(Bad_alloc("New fails "));
      MX->ReSize(10,20);
//.........这里部分代码省略.........