C++ Mat::diag方法代码示例

  double BLSSS::log_model_prob(const Selector &g)const{
    // borrowed from MLVS.cpp
    double num = vpri_->logp(g);
    if(num==BOOM::negative_infinity() || g.nvars() == 0) {
      // If num == -infinity then it is in a zero support point in the
      // prior.  If g.nvars()==0 then all coefficients are zero
      // because of the point mass.  The only entries remaining in the
      // likelihood are sums of squares of y[i] that are independent
      // of g.  They need to be omitted here because they are omitted
      // in the non-empty case below.
      return num;
    }
    SpdMatrix ivar = g.select(pri_->siginv());
    num += .5*ivar.logdet();
    if(num == BOOM::negative_infinity()) return num;

    Vector mu = g.select(pri_->mu());
    Vector ivar_mu = ivar * mu;
    num -= .5*mu.dot(ivar_mu);

    bool ok=true;
    ivar += g.select(suf().xtx());
    Mat L = ivar.chol(ok);
    if(!ok)  return BOOM::negative_infinity();
    double denom = sum(log(L.diag()));  // = .5 log |ivar|
    Vec S = g.select(suf().xty()) + ivar_mu;
    Lsolve_inplace(L,S);
    denom-= .5*S.normsq();  // S.normsq =  beta_tilde ^T V_tilde beta_tilde
    return num-denom;
  }

void train_test(int nclasses, const Mat &train_data, const Mat &train_labels, const Mat &test_data, const Mat &test_labels, Mat &confusion) {
    // setup the ann:
    int nfeatures = train_data.cols;
    Ptr<ml::ANN_MLP> ann = ml::ANN_MLP::create();
    Mat_<int> layers(4,1);
    layers(0) = nfeatures;     // input
    layers(1) = nclasses * 8;  // hidden
    layers(2) = nclasses * 4;  // hidden
    layers(3) = nclasses;      // output, 1 pin per class.
    ann->setLayerSizes(layers);
    ann->setActivationFunction(ml::ANN_MLP::SIGMOID_SYM,0,0);
    ann->setTermCriteria(TermCriteria(TermCriteria::MAX_ITER+TermCriteria::EPS, 300, 0.0001));
    ann->setTrainMethod(ml::ANN_MLP::BACKPROP, 0.0001);

    // ann requires "one-hot" encoding of class labels:
    Mat train_classes = Mat::zeros(train_data.rows, nclasses, CV_32FC1);
    for(int i=0; i<train_classes.rows; i++)
    {
        train_classes.at<float>(i, train_labels.at<int>(i)) = 1.f;
    }
    cerr << train_data.size() << " " << train_classes.size() << endl;

    ann->train(train_data, ml::ROW_SAMPLE, train_classes);

    // run tests on validation set:
    for(int i=0; i<test_data.rows; i++) {
        int pred  = ann->predict(test_data.row(i), noArray());
        int truth = test_labels.at<int>(i);
        confusion.at<int>(pred, truth) ++;
    }
    Mat correct = confusion.diag();
    float accuracy = sum(correct)[0] / sum(confusion)[0];
    cerr << "accuracy: " << accuracy << endl;
    cerr << "confusion:\n" << confusion << endl;
}

      void expectedMatDiag() {
        if(_colInd >= _genMat.n_cols) {
          return;
        }

        cout << "- Compute expectedMatDiag() ... ";
        save<double>("Mat.diagSuper", _genMat.diag(_colInd));
        cout << "done." << endl;
      }

      void expectedMatDiagElemDivide() {
        if(_rowInd >= _genMat.n_rows) {
          return;
        }

        if (_genRowVec.n_elem != min(_genMat.n_rows - _rowInd, _genMat.n_cols)) {
          return;
        }

        cout << "- Compute expectedMatDiagElemDivide() ... ";

        _genMat.diag(-_rowInd) /= _genRowVec;
        save<double>("Mat.diagSubElemDivide", _genMat);

        cout << "done." << endl;
      }

  double BM::Loglike(const Vector &ab, Vec &g, Mat &h, uint nd) const{
    if (ab.size() != 2) {
      report_error("Wrong size argument.");
    }
    double alpha = ab[0];
    double beta = ab[1];
    if (alpha <= 0 || beta <= 0) {
      if (nd > 0) {
        g[0] = (alpha <= 0) ? 1.0 : 0.0;
        g[1] = (beta <= 0) ? 1.0 : 0.0;
        if (nd > 1) {
          h = 0.0;
          h.diag() = -1.0;
        }
      }
      return negative_infinity();
    }

    double n = suf()->n();
    double sumlog = suf()->sumlog();
    double sumlogc = suf()->sumlogc();

    double ans = n*(lgamma(alpha + beta) - lgamma(alpha)-lgamma(beta));
    ans += (alpha-1)*sumlog + (beta-1)*sumlogc;

    if(nd>0){
      double psisum = digamma(alpha + beta);
      g[0] = n*(psisum-digamma(alpha)) + sumlog;
      g[1] = n*(psisum-digamma(beta)) + sumlogc;

      if(nd>1){
 	double trisum = trigamma(alpha+beta);
 	h(0,0) = n*(trisum - trigamma(alpha));
 	h(0,1) = h(1,0) = n*trisum;
 	h(1,1) = n*(trisum - trigamma(beta));}}
    return ans;
  }

  //----------------------------------------------------------------------
  double LSB::log_model_prob(const Selector &g)const{
    double num = vs_->logp(g);
    if(num==BOOM::negative_infinity()) return num;

    Ominv = g.select(pri_->siginv());
    num += .5*Ominv.logdet();
    if(num == BOOM::negative_infinity()) return num;

    Vec mu = g.select(pri_->mu());
    Vec Ominv_mu = Ominv * mu;
    num -= .5*mu.dot(Ominv_mu);

    bool ok=true;
    iV_tilde_ = Ominv + g.select(suf()->xtx());
    Mat L = iV_tilde_.chol(ok);
    if(!ok)  return BOOM::negative_infinity();
    double denom = sum(log(L.diag()));  // = .5 log |Ominv|

    Vec S = g.select(suf()->xty()) + Ominv_mu;
    Lsolve_inplace(L,S);
    denom-= .5*S.normsq();  // S.normsq =  beta_tilde ^T V_tilde beta_tilde

    return num-denom;
  }

inline
bool
op_logmat_cx::apply_common(Mat< std::complex<T> >& out, Mat< std::complex<T> >& S, const uword n_iters)
  {
  arma_extra_debug_sigprint();
  
  typedef typename std::complex<T> eT;
  
  Mat<eT> U;
  
  const bool schur_ok = auxlib::schur(U,S);
  
  if(schur_ok == false)  { arma_extra_debug_print("logmat(): schur decomposition failed"); return false; }
  
//double theta[] = { 1.10e-5, 1.82e-3, 1.62e-2,               5.39e-2,               1.14e-1,               1.87e-1,               2.64e-1              };
  double theta[] = { 0.0,     0.0,     1.6206284795015624e-2, 5.3873532631381171e-2, 1.1352802267628681e-1, 1.8662860613541288e-1, 2.642960831111435e-1 };
  // theta[0] and theta[1] not really used
  
  const uword N = S.n_rows;
  
  uword p = 0;
  uword m = 6;
  
  uword iter = 0;
  
  while(iter < n_iters)
    {
    const T tau = norm( (S - eye< Mat<eT> >(N,N)), 1 );
    
    if(tau <= theta[6])
      {
      p++;
      
      uword j1 = 0;
      uword j2 = 0;
      
      for(uword i=2; i<=6; ++i)  { if( tau      <= theta[i])  { j1 = i; break; } }
      for(uword i=2; i<=6; ++i)  { if((tau/2.0) <= theta[i])  { j2 = i; break; } }
      
      // sanity check, for development purposes only
      arma_debug_check( (j2 > j1), "internal error: op_logmat::apply_direct(): j2 > j1" );
      
      if( ((j1 - j2) <= 1) || (p == 2) )  { m = j1; break; }
      }
    
    const bool sqrtmat_ok = op_sqrtmat_cx::apply_direct(S,S);
    
    if(sqrtmat_ok == false)  { arma_extra_debug_print("logmat(): sqrtmat() failed"); return false; }
    
    iter++;
    }
  
  if(iter >= n_iters)  { arma_debug_warn("logmat(): reached max iterations without full convergence"); }
  
  S.diag() -= eT(1);
  
  if(m >= 1)
    {
    const bool helper_ok = op_logmat_cx::helper(S,m);
    
    if(helper_ok == false)  { return false; }
    }
  
  out = U * S * U.t();
  
  out *= eT(eop_aux::pow(double(2), double(iter)));
  
  return true;
  }

//.........这里部分代码省略.........
    cout << " typename(weird)=" << typeid(weird).name() << endl;
    cout << " typename(weird.real)=" << typeid(weird.real()).name() << endl;
    cout << " typename(weird.imag)=" << typeid(weird.imag()).name() << endl;

    cout << " weird.real()=" << weird.real();
    cout << " weird.imag()=" << weird.imag();

    cout << "sizeof(sym<3,cplx>)=" << sizeof(sym) << " sizeof(mat<2,3,sym>=" << sizeof(weird) << endl;

    Mat<2,3> m23(twoXthree);
    Mat<3,1> m31(threeXone);
    cout << "m23=" << m23 << endl;
    cout << "m31=" << m31 << endl;
    cout << "m23*-m31=" << m23*-m31 << endl;
    cout << "~ ~m31 * ~-m23=" << ~((~m31)*(~-m23)) << endl;

    Mat<2,3,Complex> c23(m23);
    Mat<3,1,Complex> c31(m31);
    cout << "c23=" << c23 << endl;
    cout << "c31=" << c31 << endl;
    cout << "c23*c31=" << c23*-c31 << endl;
    cout << "  ~c31 * ~-c23=" << (~c31)*(~-c23) << endl;
    cout << "~ ~-c31 * ~c23=" << ~((~-c31)*(~c23)) << endl;


    Mat<3,4> m34;
    Mat<3,4,Complex> cm34;


    cm34 = mdc;
    m34 = cm34.real();

    cout << "Mat<3,4,Complex> cm34=" << cm34 << endl;
    cout << "cm34.diag()=" << cm34.diag() << endl;

    cout << "cm34 + cm34=" << cm34+cm34 << endl; //INTERNAL COMPILER ERROR IN Release MODE
    cout << "~cm34 * 1000=" << ~cm34 * 1000. << endl;

    cout << "m34=" << m34 << endl;
    m34 =19.123;
    cout << "after m34=19.123, m34=" << m34 << endl;
 
    const double ddd[] = { 11, 12, 13, 14, 15, 16 }; 
    const complex<float> ccc[] = {  complex<float>(1.,2.),  
                                    complex<float>(3.,4.),
                                    complex<float>(5.,6.),
                                    complex<float>(7.,8.) };
    Vec<2,complex<float>,1> cv2(ccc);
    cout << "cv2 from array=" << cv2 << endl;
    cv2 = Vec<2,complex<float> >(complex<float>(1.,2.), complex<float>(3.,4.));
    cout << "cv2 after assignment=" << cv2 << endl;

    cout << "cv2.real()=" << cv2.real() << " cv2.imag()=" << cv2.imag() << endl;

    Vec<2,negator<complex<float> >,1>& negCv2 = (Vec<2,negator<complex<float> >,1>&)cv2;
    Vec<2,conjugate<float>,1>& conjCv2 = (Vec<2,conjugate<float>,1>&)cv2;
    Vec<2,negator<conjugate<float> >,1>& negConjCv2 = (Vec<2,negator<conjugate<float> >,1>&)cv2;

    

    Vec<2,complex<float> > testMe = cv2;
    cout << "testMe=cv2 (init)=" << testMe << endl;
    testMe = cv2;
    cout << "testMe=cv2 (assign)=" << testMe << endl;

    //Class constructor
	TimeSegmentation(Tobj &G, Col <T1> map_in, Col <T1> timeVec_in, uword a, uword b, uword c, uword interptype = 1,
					 uword shots = 1)
    {
		cout << "Entering Class constructor" << endl;
	    n1 = a; //Data size
	    n2 = b;//Image size
	    L = c; //number of time segments
	    type = interptype; // type of time segmentation performed
	    Nshots = shots; // number of shots
	    obj = &G;
	    fieldMap = map_in;
    	cout << "N1 = " << n1 << endl;
		cout << "N2 = " << n2 << endl;
		cout << "L = " << L << endl;


		AA.set_size(n1, L); //time segments weights
	    timeVec = timeVec_in;
	    T_min =timeVec.min();
	    T1 rangt = timeVec.max()-T_min;
		tau = (rangt + datum::eps) / (L - 1); // it was L-1 before
	    timeVec = timeVec-T_min;

	    uword NOneShot = n1/Nshots;
	    if (L==1) {
		    tau = 0;
		    AA.ones();
	    }
	    else {
			Mat <CxT1> tempAA(NOneShot, L);
		    if (type==1) {// Hanning interpolator
			    cout << "Hanning interpolation" << endl;
			    for (unsigned int ii = 0; ii<L; ii++) {
				    for (unsigned int jj = 0; jj<NOneShot; jj++) {
					    if ((std::abs(timeVec(jj)-((ii)*tau)))<=tau) {
						    tempAA(jj, ii) = 0.5+0.5*std::cos((datum::pi)*(timeVec(jj)-((ii)*tau))/tau);
					    }
					    else {
						    tempAA(jj, ii) = 0.0;
					    }
				    }
			    }
			    AA = repmat(tempAA, Nshots, 1);
		    }
		    else if (type==2) { // Min-max interpolator: Exact LS interpolator

			    cout << "Min Max time segmentation" << endl;

			    Mat <CxT1> Ltp;
			    Ltp.ones(1, L);
			    Col <CxT1> ggtp;
			    ggtp.ones(n2, 1);
			    Mat <CxT1> gg;
			    gg = exp(i*fieldMap*tau)*Ltp;
			    Mat <CxT1> iGTGGT;
			    iGTGGT.set_size(L+1, n2);
			    Mat <CxT1> gl;
			    gl.zeros(n2, L);


			    for (unsigned int ii = 0; ii<L; ii++) {
				    for (unsigned int jj = 0; jj<n2; jj++) {
					    gl(jj, ii) = pow(gg(jj, ii), (T1) (ii+1));
				    }
			    }

			    Mat <CxT1> G;
			    G.set_size(n2, L);

			    for (unsigned int jj = 0; jj<L; jj++) {
				    if (jj==0) {
					    G.col(jj) = ggtp;
				    }
				    else {
					    G.col(jj) = gl.col(jj-1);
				    }
			    }

			    Col <CxT1> glsum;
			    Mat <CxT1> GTG;
			    GTG.zeros(L, L);
			    GTG.diag(0) += n2;
			    glsum = sum(gl.t(), 1);
				Mat <CxT1> GTGtp(L, L);
				for (unsigned int ii = 0; ii < (L - 1); ii++) {
					GTGtp.zeros();
				    GTGtp.diag(-(T1) (ii+1)) += glsum(ii);
				    GTGtp.diag((T1) (ii+1)) += std::conj(glsum(ii));
				    GTG = GTG+GTGtp;
			    }

			    T1 rcn = 1/cond(GTG);
			    if (rcn>10*2e-16) { //condition number of GTG
				    iGTGGT = inv(GTG)*G.t();

			    }
			    else {
				    iGTGGT = pinv(GTG)*G.t(); // pseudo inverse
			    }
//.........这里部分代码省略.........