C++ mat::insert_cols方法代码示例

mat addPadding(mat x, int ksize) {
    int offset = ksize/2;
    x.insert_rows(0, offset);
    x.insert_rows(x.n_rows, offset);
    x.insert_cols(0, offset);
    x.insert_cols(x.n_cols, offset);
    return(x);
}

void SolveX(mat& C,mat& D,mat& G,mat& B,mat& A,mat& X,vec& Z){//Solve vector of unknowns
    A.zeros(0,0);
    C = trans(B);
    A.insert_rows(0,G);
    A.insert_cols(A.n_cols,B);
    C.insert_cols(C.n_cols,D);
    A.insert_rows(A.n_rows,C);
    X = solve( A, Z );
}

bool inverseMatrix(const mat& A, mat& inv)
{
	unsigned int n = A.n_rows;
	mat LU(n, n);
	vector<unsigned int> Pvec;
	if (LUP(A, LU, Pvec) == false)
		return false;

	for (unsigned int i = 0; i < n; ++i)
	{
		vec b = zeros<vec>(n);
		b(distance(Pvec.begin(), find(Pvec.begin(), Pvec.end(), i))) = 1;
		inv.insert_cols(i, SLAU(LU, b));
	}
	return true;
}

//.........这里部分代码省略.........
            // construct normalization constant, q0i, to sample s(i)
            // build loqq0 and exponentiate
            colvec bki(T), bbar_ki(T); /* T x 1, D_k^-1*Omega_k*b_ki = C(k,0)*b_ki */
            mat bbar_i(K,T); bbar_i.zeros();
            double logd_dk = 0; /* set of T 0 mean gaussian densities for term k */
            double logq0ki = 0, logq0i = 0, q0i = 0;
            // accumulate weight, q0i, for s(i) over K iGMRF terms  
            for( k = 0; k < K; k++)
            {
                 logq0ki       = 0; /* reset k-indexed log-like on each k */
                 //a1k           = 0.5*(double(T)) + a;
                 a1k           = 0.5*(double(T)-double(o_adjust(k))) + a;
                 bki           = B.slice(k).row(i).t();
                 bbar_ki       = C(k,0) * bki; /* T x 1 */
                 bbar_i.row(k) = bbar_ki.t();
                 B1(k,i)       = 0.5*dot( D.row(k), pow((bki-bbar_ki),2) ); /* no b */
                 logd_dk       = 0; /* set of T gaussian densities for term k */
                 /* dmvn(zro|m,Q.slice(k),true) */
                 for( j = 0; j < T; j++ )
                 {
                    logd_dk   += R::dnorm(0.0,0.0,double(1/sqrt(D(k,j))),true);
                 }
                 logq0ki      = logd_dk + lgamma(a1k) + a*log(b) -
                                   lgamma(a) - a1k*trunc_log(B1(k,i)+b);
                 logq0i       += logq0ki;
            } /* end loop k over iGMRF terms */
            q0i = trunc_exp(logq0i);

            // construct posterior sampling weights to sample s(i)
            colvec weights(M+1); weights.zeros();
            /* evaluate likelihood under kappa_star(k,i) */
            double lweights_l;
            for(l = 0; l < M; l++) /* cycle through all clusters for s(i) */
            {
                s(i)          = l; /* will compute likelihoods for every cluster */  
                lweights_l = 0; /* hold log densities for K computations */
                for(k = 0; k < K; k++)
                {
                    bki            = B.slice(k).row(i).t();
                    for( j = 0; j < T; j++ )
                    {
                      /* effectively making assignment, s(i) = l */
                      lweights_l   += trunc_log(R::dnorm(bki(j),bbar_i(k,j),
                                    double(1/sqrt(kappa_star(k,l)*D(k,j))),false));
                    } /* end loop j over time index */
                } /* end loop k over iGMRF terms */
                //if(lweights_l < -300){lweights_l = -300;}
                weights(l)          = trunc_exp(lweights_l);
                weights(l)          *= double(Num(s(i)))/(double(N) - 1/ipr(i) + conc);
            } /* end loop l over existing or populated clusters */
            /* M+1 or new component sampled from F_{0} */
            weights(M)              = conc/(double(N) - 1/ipr(i) + conc)*q0i;

            // normalize weights
            sweights = sum(weights);
            if(sweights == 0)
            {
                weights.ones(); weights *= 1/(double(M)+1);
            }
            else
            {
                weights /= sweights;
            }

            // conduct discrete posterior draw for s(j)
            unsigned long MplusOne = M + 1;
            s(i) = rdrawone(weights, MplusOne);

            // if new cluster chosen, generate new location
            if(s(i) == M)
            {
                /* sample posterior of ksi_star[k,m] for 1 (vs. n_m) observation */
                double a_star_k; /* shape for 1 obs */
                double bstar_ki;
                kappa_star.insert_cols(M,1); /* add K vector new location to kappa_star */
                num.insert_rows(M,1);
                Num.insert_rows(M,1);
                for(k = 0; k < K; k++)
                {
                     a_star_k         = 0.5*(double(T) - double(o_adjust(k))) + a; /* shape for 1 obs */
                     bstar_ki         = B1(k,i) + b; /* B1(k,i) is a scalar quadratic product */
                     /*
                     bki              = B.slice(k).row(i).t();
                     bstar_ki         = 0.5*( as_scalar(bki.t()*symmatl(Q.slice(k))*bki) ) + b;
                     */
                     kappa_star(k,M)  = rgamma(1, a_star_k, (1/bstar_ki))[0];
                }
                num(M)   = 1;
                Num(M)   = 1/ipr(i);
                M        = MplusOne;
            }
            else
            {
                num(s(i)) += 1;
                Num(s(i)) += 1/ipr(i);
            }
            
        } /* end loop i for cluster assignment to unit i = 1,...,N */
        END_RCPP
    } /* end function bstep for cluster assignments, s, and computing zb */