C++ cube::slice方法代码示例

cube barrel_distort_rgb(const cube &F, double offset_x) {
  cube G(F.n_rows, F.n_cols, F.n_slices);
  G.slice(0) = barrel_distort(F.slice(0), offset_x);
  G.slice(1) = barrel_distort(F.slice(1), offset_x);
  G.slice(2) = barrel_distort(F.slice(2), offset_x);
  return G;
}

SEXP move_B(const mat& y, cube& B, const mat& kappa_star, const field<mat>& C,
            //const field<sp_mat>& C,
            mat& gamma, const mat& D, const ucolvec& s, double tau_e)
{
     BEGIN_RCPP
     // N x T matrix of standardized counts, y
     // B = (B_1,...,B_K), where B_k is N x T matrix for iGMRF term k
     // N x T matrix, gamma = sum_{k=1^K}(B_k)
     // K x T, D, where row k contains T x 1 diagonal elements of Q_k
     // sample cluster assignments, s(1), ..., s(N)  
     // Q = (Q_1,...,Q_K), where Q_k is a T x T de-scaled iGMRF precision matrix
     // C = (C_1,...,C_K), where C_k = D_k^-1 * Omega_k, 
     // where Omega_k is the T x T adjacency matrix for iGMRF term, k
     // D is a K x T matrix where row k contains T diagonal elements of Q_k
     // K x M matrix, kappa_star records locations for each iGMRF term 
     int K     = B.n_slices;
     int N     = y.n_rows;
     int T     = D.n_cols;
     colvec bbar_ki(T); bbar_ki.zeros();
     rowvec gammatilde_ki(T); gammatilde_ki.zeros();
     rowvec ytilde_ki(T); ytilde_ki.zeros();
     rowvec d_k(T);
     vec zro(T); zro.zeros();
     double e_ij, phi_ij, bkij, h_ij;
     int k = 0, i = 0, j = 0; /* loop variables */
     for( k = 0; k < K; k++ ) /* over iGMRF terms */
     {
          d_k                           = D.row(k);
          for( i = 0; i < N; i++ ) /* over units */
          {
               /* take out T x 1, b_ki, to be sampled from gamma */
               //gammatilde_ki        = gamma.row(i) - B.slice(k).row(i); /* 1 x T */
               //ytilde_ki            = y.row(i) - gammatilde_ki;
               gammatilde_ki            = gamma.row(i); /* 1 x T */
               // mean of univariate iGMRF, b_kij = 1/d_kj * (omega_kj(-j) * b_ki(-j))
               bbar_ki                  = C(k,0) * B.slice(k).row(i).t(); /* T x 1 */
               for( j = 0; j < T; j++ ) /* over time points */
               {
                    gammatilde_ki(j)    -= B.slice(k)(i,j);
                    ytilde_ki(j)        = y(i,j) - gammatilde_ki(j);
                    // mean of univariate iGMRF, b_kij = 1/d_kj * (omega_kj(-j) * b_ki(-j))
                    B.slice(k)(i,j)     = 0;
                    //bbar_ki(j)          = dot( C(k,0).row(j), B.slice(k).row(i) );
                    e_ij                = tau_e*ytilde_ki(j) 
                                             + d_k(j)*kappa_star(k,s(i)) * bbar_ki(j);
                    phi_ij              = tau_e + d_k(j)*kappa_star(k,s(i));
                    h_ij                = (e_ij / phi_ij);
                    bkij                = rnorm( 1, (h_ij), sqrt(1/phi_ij) )[0];
                    B.slice(k)(i,j)     = bkij;
                    /* put back new sampled values for b_kij */
                    gammatilde_ki(j)    += B.slice(k)(i,j);
               } /* end loop j over time points */
               /* put back new sampled values for b_ki */
               //gamma.row(i)         = gammatilde_ki + B.slice(k).row(i); 
               gamma.row(i)        = gammatilde_ki;
          } /* end loop i over units */
     } /* end loop k over iGMRF terms */
     END_RCPP
} /* end function ustep to sample B_1,...,B_K */

/**********************************************************************
*   MSE Class
***********************************************************************/
double MSE::cost(cube pred, cube y){
    dbg_assert(y.n_cols == 1 && y.n_slices == 1);
    pred.reshape(pred.n_elem,1,1);
    y.reshape(y.n_elem,1,1);
    
    double retn = 0.0f;
    for(int i=0; i < y.n_elem; i++) {
        retn += pow(pred(i,0,0) - y(i,0,0),2);
    }
    retn /= y.n_elem;
    retn *= 0.5f;
    //return(retn);
    return(0.5 * mean(mean(square(pred.slice(0) - y.slice(0)))));        
}

cube Utils::conv3(cube body, cube kernel)
{
	if (body.max() <= 0)
	{
		return zeros(body.n_rows, body.n_cols, body.n_slices);
	}
	cube Q(body);
	for (int i = 1; i < body.n_slices - 1; i++)
	{
		Q.slice(i) = conv2(body.slice(i), kernel.slice(1), "same") +
			conv2(body.slice(i - 1), kernel.slice(0), "same") +
			conv2(body.slice(i + 1), kernel.slice(2), "same");
	}
	return Q;
}

cube noncont_slices(const cube & X, const uvec & index){
	cube Xsub(X.n_rows, X.n_cols, index.n_elem);
	for(unsigned int i=0; i<index.n_elem; i++){
		Xsub.slice(i) = X.slice(index[i]);
	}
	return Xsub; 
}

void lmMean(const cube& wX, const rowvec& vec_beta, mat& mean){
	int n = wX.n_slices, P=wX.n_cols;
	mean = zeros(n,P);
	for(int i=0;i<n;i++){
		mean.row(i) = vec_beta* wX.slice(i);
	}
}

cube conv2(const cube &F, const mat &H) {
  cube G(F.n_rows, F.n_cols, F.n_slices);
  for (uword i = 0; i < F.n_slices; i++) {
    G.slice(i) = conv2(F.slice(i), H);
  }
  return G;
}

cube imresize2(const cube &C, uword m, uword n) {
  cube F(m, n, C.n_slices);
  for (uword k = 0; k < C.n_slices; k++) {
    F.slice(k) = imresize2(C.slice(k), m, n);
  }
  return F;
}

//! compute (trial)point-to-(model)point Mahalanobis distance
//! @param[in] index		index of the points (in trial and model) to be compared
//! @param[in] &trial		reference to the trial
//! @param[in] &model		reference to the model
//! @param[in] &variance	reference to the model variance
//! @return 			    Mahalanobis distance between trial-point and model-point
float Classifier::mahalanobisDist(int index,mat &trial,mat &model,cube &variance)
{
	mat difference = trial.col(index) - model.col(index);
	mat distance = (difference.t() * (variance.slice(index)).i()) * difference;

	return distance(0,0);
}

void white(const mat& y, const mat& y_center, const uvec& sti, const cube& matphi1, const cube& matphi2, mat& wy, mat& wy_center, const uvec& WTIME){
	int L=matphi1.n_slices, wT = WTIME.n_elem, t=0; 
	mat multiply;
	if(L==0){wy=y; wy_center = y_center;}
	else{
		wy = y.rows(WTIME);
		wy_center = y_center.rows(WTIME);
		for(int wt=0; wt<wT ; wt++){
			t = WTIME[wt];
			for(int l=0; l<L; l++){
				multiply = (sti[t-l-1])?(matphi1.slice(l)):(matphi2.slice(l));
				wy.row(wt) -= y.row(t-l-1)*multiply;					
				wy_center.row(wt) -= y_center.row(t-l-1)*multiply;
			}
		}
	}
}

cube ConvLayer::backward(cube delta) {
   
    // NOTE: delta may come from a linear layer
    delta.reshape(_a.n_rows, _a.n_cols, _a.n_slices);
    
    if(!_prev) {
        dbg_print("null pointer to _prev in ConvLayer::backward");
        return zeros<cube>(0,0,0);
    }

    // Compute Weight Updates
    cube input = addPadding(_prev->getActivationCube(), _ksize);
    for(int i = 0; i < _units; i++){
        // dz_dw
        for(int j=0; j < input.n_slices; j++) {
            _dw(i).slice(j) = conv2d(input.slice(j), delta.slice(i));
        }
        // dz_db
        _db(0,0,i) = accu(delta.slice(i));
    }
    

    // Compute next delta
    int nr = _prev->getActivationCube().n_rows;
    int nc = _prev->getActivationCube().n_cols;
    int ns = _prev->getActivationCube().n_slices;
    
    delta = addPadding(delta, _ksize);
    cube next_delta = zeros<cube>(nr,nc,ns);
    for(int i=0; i < ns; i++){
        for(int j=0; j < _units; j++) {
          next_delta.slice(i) += conv2d(
                                    delta.slice(j), 
                                    fliplr(flipud(_w(j).slice(i)))
                                    );
        }
    }

    if(_prev) {
        return(_prev->backward(next_delta));
    }

    return(next_delta);

}

void lm2(const mat& wy, const cube& wX, const mat& Omega, rowvec& numer, mat& denom){
	int P=wy.n_cols, n = wy.n_rows, k = wX.n_rows;
	denom.zeros(), numer.zeros();
	mat twX(P,k);
	for(int t=0; t<n; t++){
		twX = trans(wX.slice(t));
		denom=denom+wX.slice(t)*Omega*twX; 
		numer=numer+wy.row(t)*Omega*twX;
	}
}

cube maxPoolLayer::forward (cube x) {

    for(int i=0; i < x.n_slices; i++)
        _a.slice(i) = maxDownSample(x.slice(i), 
                                    _sample_size, 
                                    _mask.slice(i));
    
    if(_next)
        return(_next->forward(_a));
    return _a;
}

SEXP move_kappastar_alt(mat& kappa_star, const cube& B, const cube& Q, 
                    const ucolvec& s, uvec& o, 
                    int T, int a, int b, const vec& ipr)
{
     BEGIN_RCPP
     // K x M matrix, kappa_star, records location values 
     // N x T matrices, B_1,...,B_K contains the de-noised functions
     // on posterior for kappa_star.  
     int K = kappa_star.n_rows; /* number of iGMRF terms */
     int M = kappa_star.n_cols; /* number of clusters */
     int k = 0, m = 0, i = 0, count_m;
     double num_m;
     uvec pos_m;
     double a1_mk; /* posterior shape parameter */
     double b1_mk; /* posterior rate parameter */
     colvec b_ki(T);
     for(m = 0; m < M; m++)
     {
          pos_m          = find( s == m ); /* s is vector length N */
          count_m        = pos_m.n_elem;
          num_m          = sum( 1/ipr(pos_m) );
          /* sample posterior for kappa_star(k,) for each iGMRF term, k = 1,...,K */
          for( k = 0; k < K; k++ )
          {
               b1_mk           = 0;
               for( i = 0; i < count_m; i++ )
               {
                    b_ki                  = B.slice(k).row(pos_m(i)).t(); /* T x 1 */
                    b1_mk                 += 0.5*( as_scalar(b_ki.t()*symmatl(Q.slice(k))*b_ki)
                                                       / ipr(pos_m(i)) );       
               } /* end loop i over num_m weighted units in cluster m */
               b1_mk               += b;
               a1_mk               = 0.5*num_m*(double(T)-double(o(k))) + a;
               kappa_star(k,m)     = rgamma(1, a1_mk, (1/b1_mk))[0];
          } /* end loop over K iGMRF terms */ 
          
     } /* end loop m over clusters */
     /* add a bumper */
     END_RCPP
} /* end function move_kappastar() to sample cluster locations */

void phi_design(mat& U, const uvec& sti, const int L, umat& vxicat, cube& wZ, const uvec& WTIME){
	int P=U.n_cols, wT = WTIME.n_elem, t;
	int Psq=P*P;
	wZ.zeros();
	for(int wt=0; wt<wT; wt++){
		t = WTIME[wt];
		for(int l=1;l<=L;l++){
			if(sti[t-l]){
				for(int p=0;p<P;p++){
					uvec cond= (vxicat.col(p) >= l);
					wZ.slice(wt).col(p).subvec((l-1)*Psq+ p*P,(l-1)*Psq+(p+1)*P-1) = (cond)%trans(U.row(t-l));
				}
				wZ.slice(wt).rows((L+l-1)*Psq,(L+l)*Psq-1).zeros();
			}else{
				wZ.slice(wt).rows((l-1)*Psq,l*Psq-1).zeros();
				for(int p=0;p<P;p++){
					uvec cond= (vxicat.col(p) >= l);
					wZ.slice(wt).col(p).subvec((L+l-1)*Psq+p*P,(L+l-1)*Psq+(p+1)*P-1)=cond%trans(U.row(t-l));
				}
			}
		}
	}

}