C++ NumericMatrix::row方法代码示例

// [[Rcpp::export]]
NumericVector loopC(NumericVector nodelist,int al,IntegerVector x, IntegerVector x2, List pm_lc, NumericMatrix Lix, int finind){
  NumericVector::iterator k;
  LogicalVector res;
  IntegerVector daughter;
  int temp1, temp2;
  int k2;
  int l = 0;
  int g = 0;
  int n = x.size();
  LogicalVector in1(x2.size());
  LogicalVector in2(x2.size());
  IntegerVector indices = seq_len(n); //from 1 to n...
  for(k = nodelist.begin(); k != nodelist.end(); ++k) {
    k2  = nodelist[l];
    res = x == k2;
    daughter = x2[res];
    in1 = x2==daughter[0];
    in2 = x2==daughter[1];
    temp1 = as<int>(indices[in1]);
    temp2 = as<int>(indices[in2]);
    NumericMatrix pmtmp1 = pm_lc[temp1 - 1];
    NumericMatrix pmtmp2 = pm_lc[temp2 - 1];
    for(g = 0; g<al; ++g){
      Lix(k2 - 1, g) = sum(Lix.row(daughter[0] - 1) * pmtmp1.row(g)) *
      sum(Lix.row(daughter[1] - 1) * pmtmp2.row(g));
    }
    l=l+1;
  }
  NumericVector res2 = Lix(finind,_);
  return res2;
}

// [[Rcpp::export]]
List runit_Row_Column_sugar( NumericMatrix x){
    NumericVector r0 = x.row(0) ;
    NumericVector c0 = x.column(0) ;
    return List::create(
        r0,
        c0,
        x.row(1),
        x.column(1),
        x.row(1) + x.column(1)
        ) ;
}

// [[Rcpp::export]]
NumericMatrix max_ll_matrix(NumericMatrix x, NumericMatrix y, IntegerVector which_max) {
  NumericMatrix out(x.nrow(), x.ncol());
  for (int i=0; i < x.nrow(); i++) {
    if (which_max[i] == 0)
      out.row(i) = x.row(i);
    else if (which_max[i] == 1)
      out.row(i) = y.row(i);
    else
      std::cout << "ERROR: which_max values must be in (0, 1)" << std::endl;
  }
  return out;
}

// [[Rcpp::export]]
NumericMatrix imp_neighbour_avg(NumericMatrix x, double k) {
  // input matrix is expected to have >= 3 columns
  NumericMatrix ans = clone(x);
  int nr = ans.nrow(), nc = ans.ncol();

  for(int i = 0; i < nr; i++) {
    // first and last values are set to 0 if NA
    if (R_IsNA(ans(i, 0))) ans(i, 0) = k;
    if (R_IsNA(ans(i, nc-1))) ans(i, nc-1) = k;
    
    for(int j = 1; j < (nc-1); j++) {
      if (R_IsNA(ans(i,j))) {
	// if the next value is NA and all previous values are 0
	// then we set to 0	
	if (R_IsNA(ans(i,j+1))) {
	  NumericVector v = subset(ans.row(i), j);
	  if (allZero(v, k)) ans(i,j) = k;
	} else { // next is not NA, set to mean of neighbours
	  ans(i,j) = (ans(i,j-1) + ans(i,j+1))/2;
	}
      }
    }
  }    
  return(ans);
}

// [[Rcpp::export]]
NumericMatrix rcpp_js_distance(NumericMatrix mat) {
  
   // allocate the matrix we will return
   NumericMatrix rmat(mat.nrow(), mat.nrow());
   
   for (int i = 0; i < rmat.nrow(); i++) {
      for (int j = 0; j < i; j++) {
      
         // rows we will operate on
         NumericMatrix::Row row1 = mat.row(i);
         NumericMatrix::Row row2 = mat.row(j);
         
         // compute the average using std::tranform from the STL
         std::vector<double> avg(row1.size());
         std::transform(row1.begin(), row1.end(), // input range 1
                        row2.begin(),             // input range 2
                        avg.begin(),              // output range 
                        average);                 // function to apply
      
         // calculate divergences
         double d1 = kl_divergence(row1.begin(), row1.end(), avg.begin());
         double d2 = kl_divergence(row2.begin(), row2.end(), avg.begin());
        
         // write to output matrix
         rmat(i,j) = std::sqrt((double)(.5 * (d1 + d2)));
      }
   }
   
   return rmat;
}

//' @title Calculate \eqn{L_q} distance
//' @description Calculate \eqn{L_q} distance of all vectors in a matrix to a reference
//' vector.
//' @param x A numeric matrix Missing values are allowed.
//' @param ref An integer specifying the reference row.
//' @param q An integer specifying the which norm to take the L-q distance of.
//' @return a numeric vector of length \code{nrow(x) - 1}
// [[Rcpp::export]]
NumericVector dist_q_matrix (NumericVector& x_ref, NumericMatrix& x_rest, int& q) {
  int nr = x_rest.nrow();
  NumericVector out(nr);
  for (int k = 0; k < nr; k++) {
    out[k] = dist_q(x_ref, x_rest.row(k), q);
  }
  return out;
}

// [[Rcpp::export]]
NumericMatrix  genSizeComp(NumericMatrix VulnN, NumericVector CAL_binsmid, NumericMatrix selCurve,
                           double CAL_ESS, double CAL_nsamp,
                           NumericVector Linfs, NumericVector Ks, NumericVector t0s,
                           double LenCV, double truncSD) {
  int nyears = VulnN.nrow();
  int k = VulnN.ncol();
  int nbins = CAL_binsmid.size();
  NumericMatrix CAL(nyears, nbins);
  double width = CAL_binsmid(1) - CAL_binsmid(0);
  double origin = CAL_binsmid(0) - 0.5* width;
  NumericVector temp(k);
  NumericVector varAges = NumericVector::create(-0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5); // monthly ages
  for (int yr=0; yr < nyears; yr++) {
    NumericVector Nage = (VulnN.row(yr)); // numbers of catch-at-age this year
    double Ncatch = sum(Nage); // total catch this year
    if (Ncatch>0) {
      NumericVector Nage2 = (Nage/Ncatch) * CAL_ESS; // number-at-age effective sample
      List Lens(k*12);
      int count = 0;
      for (int age=1; age <= k; age++) { // loop over 1:maxage
        int Nage3 =  round(Nage2(age-1)); // number at this age
        NumericVector rands = RcppArmadillo::sample(NumericVector::create(0,1,2,3,4,5,6,7,8,9,10,11), Nage3, TRUE, NumericVector::create()) ; //  assume ages are uniformly distributed across months
        NumericVector subAgeVec = get_freq(rands, 1, 0, 12); // distribute n across months
        // NumericVector subAgeVec = Nage3/NumericVector::create(1,2,3,4,5,6,7,8,9,10,11,12); // distribute n across months
        for (int subage=0; subage<=11; subage++) { // loop over 12 months
          if (subAgeVec(subage) > 0) {
            double sage = varAges(subage) + age;
            double mean = Linfs(yr) * (1-exp(-Ks(yr)* (sage - t0s(yr)))); // calculate mean length at sub-age;
            if (mean < 0) mean = 0.01;
            NumericVector dist = tdnorm((CAL_binsmid-mean)/(LenCV*mean), -truncSD, truncSD); // prob density of lengths for this age
            NumericVector newdist = dist * selCurve(_,yr); // probability = dist * size-selection curve
            if (sum(newdist)!=0) {
              newdist = newdist/sum(newdist);
              Lens(count) = RcppArmadillo::sample(CAL_binsmid, subAgeVec(subage), TRUE, newdist); // sample lengths for this sub-age class
            } else {
              Lens(count) = NA_INTEGER;
            }
          } else {
            Lens(count) = NA_INTEGER;
          }
          count += 1;
        }
      }
      NumericVector LenVals = combine(Lens); // unlist
      NumericVector templens = get_freq(LenVals, width, origin, nbins); // calculate frequencies
      double rat = CAL_nsamp/sum(templens);
      templens =  templens * rat; // scale to CAL_nsamp
      CAL(yr,_) = templens;
    } else {
      NumericVector zeros(nbins);
      CAL(yr,_) = zeros;
    }
  }
  return(CAL);
}

// [[Rcpp::export]]
NumericVector row_kth(NumericMatrix toSort, int k) {
  int n = toSort.rows();
  NumericVector meds = NumericVector(n);
  for (int i = 0; i < n; i++) {
    NumericVector curRow = toSort.row(i);
    std::nth_element(curRow.begin(), curRow.begin() + k, curRow.end());
    meds[i] = curRow[k];
  }

  return meds;
}

// [[Rcpp::export]]
NumericMatrix knn(NumericMatrix train, IntegerVector group, NumericMatrix test, int kn, int disttype){

  int ngroup = sort_unique(group).size();

  NumericMatrix pred(test.nrow(),ngroup);
  for(int i = 0; i < test.nrow(); ++i){ // For each test point
    NumericVector dists(kn);
    dists = dists + 1000000000.0;
    IntegerVector indx(kn);

      
    for(int j = 0; j < train.nrow(); ++j){ // Run through train points
      double d = distance(train.row(j),test.row(i),disttype); // Find the distance
      if(d < dists(kn-1)){ // If distance is more than the kn furthest away (so far) it can not be one of the kn nearest neighbours. Otherwise insert in list.
	for (int k = 0; k < kn; ++k){ // Run trough list of neighbours 
	  if (d < dists(k)){// If the training point is closer than the one we are looking at, it should be inserted here
	    for (int k1 = kn-1; k1 > k; --k1){ // All neighbours further away should move one place up in the list
	      dists(k1) = dists(k1 - 1);
	      indx(k1) = indx(k1 - 1);
	    }
	    //Finally insert the new point
	    dists(k) = d;
	    indx(k) = j;
	    break;
	  }
	}
      }
    }
    //Run through neighbours
    for(int j = 0; j < indx.size(); ++j){
      // Add a vote to the group it belongs to
      pred(i,group(indx(j))) += 1.0/(double)kn;
    }
  }

  return pred;

}

NumericMatrix row_erase (NumericMatrix& x, IntegerVector& rowID) {
  rowID = rowID.sort();

  NumericMatrix x2(Dimension(x.nrow()- rowID.size(), x.ncol()));
  int iter = 0;
  int del = 1; // to count deleted elements
  for (int i = 0; i < x.nrow(); i++) {
    if (i != rowID[del - 1]) {
      x2.row(iter) = x.row(i);
      iter++;
    } else {
      del++;
    }
  }
  return x2;
}

// [[Rcpp::export]]
NumericVector row_medians(NumericMatrix toSort) {
  int n = toSort.rows();
  int medN  = toSort.cols();
  NumericVector meds = NumericVector(n);
  for (int i = 0; i < n; i++) {
    NumericVector curRow = toSort.row(i);
    std::nth_element(curRow.begin(), curRow.begin() + curRow.size()/2 - 1, curRow.end());
    double med1 = curRow[curRow.size()/2 - 1];
    if (medN % 2 == 0) {
      std::nth_element(curRow.begin(), curRow.begin() + curRow.size()/2, curRow.end());
      double med2 = curRow[curRow.size()/2];
      meds[i] = (med1 + med2)/2.0;
    } else {
      meds[i] = med1;
    }
  }

  return meds;
}

// [[Rcpp::export]]
void distance(NumericVector is, NumericVector js, NumericVector xs, NumericMatrix data, Function callback) {
  for (int i=0; i < is.length(); i++) {
    xs[i] = sqrt(sum(pow(data.row(is[i]) - data.row(js[i]), 2)));
    if (i > 0 && i % 1000 == 0) callback(1000);
  }
};

// [[Rcpp::export]]
double runit_NumericMatrix_row( NumericMatrix m){
    NumericMatrix::Row first_row = m.row(0) ;
    return std::accumulate( first_row.begin(), first_row.end(), 0.0 ) ;
}

// [[Rcpp::export]]
NumericMatrix viterbi_ths(NumericVector &theta, NumericMatrix &data,
	                  NumericVector &integrControl) {
  // theta lambda0, lambda1, lambda2, sigma, p
  // data diff of t and x
  int n = data.nrow(); int dim = data.ncol() - 1;
  double lambda0 = theta[0], lambda1 = theta[1], lambda2 = theta[2];
  double p = theta[4];
  if (lambda1 < lambda2) return NA_REAL;
  double ps0 = 1. / lambda0 / (1. / lambda0 + p / lambda1 + (1 - p) / lambda2);
  double ps1 = p / lambda1 / (1. / lambda0 + p / lambda1 + (1 - p) / lambda2);
  double ps2 = (1 - p) / lambda2 / (1. / lambda0 + p / lambda1 + (1 - p) / lambda2);
  NumericVector tt = data.column(0);
  NumericMatrix x = data(Range(0, n - 1), Range(1, dim));

  // result matrix: three cols stand for Viterbi prob of state 0,1,2
  //                at current time points. For numerical reason,
  //                the log-prob is returned.
  NumericMatrix result(n + 1, 3);
  result(0, 0) = log(ps0); result(0, 1) = log(ps1); result(0, 2) = log(ps2);
  NumericVector cartV = result.row(0);
  NumericVector cartW(3);

  // calculate all h functions
  NumericVector
    hresult00 = ths_h00(x, tt, theta, integrControl),
    hresult01 = ths_h01(x, tt, theta, integrControl),
    hresult02 = ths_h02(x, tt, theta, integrControl),
    hresult10 = ths_h10(x, tt, theta, integrControl),
    hresult11 = ths_h11(x, tt, theta, integrControl),
    hresult12 = ths_h12(x, tt, theta, integrControl),
    hresult20 = ths_h20(x, tt, theta, integrControl),
    hresult21 = ths_h21(x, tt, theta, integrControl),
    hresult22 = ths_h22(x, tt, theta, integrControl);
  
  for (int i = 0; i < n; i++) {
    NumericVector crow = x.row(i);
    if (is_true(all(crow == 0.))) {
      hresult00[i] = 0.;
      hresult01[i] = 0.;
      hresult02[i] = 0.;
      hresult10[i] = 0.;
      hresult11[i] = exp(-lambda1 * tt[i]);
      hresult12[i] = 0.;
      hresult20[i] = 0.;
      hresult21[i] = 0.;
      hresult22[i] = exp(-lambda2 * tt[i]);
    }
  }

  // calculate Viterbi path
  for (int i = 1; i <= n; i++) {
    cartW[0] = cartV[0] + log(hresult00[i - 1]);
    cartW[1] = cartV[1] + log(hresult10[i - 1]);
    cartW[2] = cartV[2] + log(hresult20[i - 1]);
    result(i, 0) = max(cartW);

    cartW[0] = cartV[0] + log(hresult01[i - 1]);
    cartW[1] = cartV[1] + log(hresult11[i - 1]);
    cartW[2] = cartV[2] + log(hresult21[i - 1]);
    result(i, 1) = max(cartW);

    cartW[0] = cartV[0] + log(hresult02[i - 1]);
    cartW[1] = cartV[1] + log(hresult12[i - 1]);
    cartW[2] = cartV[2] + log(hresult22[i - 1]);
    result(i, 2) = max(cartW);

    cartV = result.row(i);
  }

  return result;
}

// [[Rcpp::export]]
NumericMatrix partial_viterbi_ths(NumericVector &theta, NumericMatrix &data,
				  NumericVector &integrControl,
				  int &startpoint, int &pathlength){
  // theta lambda0, lambda1, lambda2, sigma, p
  // data diff of t and x
  // startpoint the start time point, note that
  //            the first time point in data is t0
  // pathlength the length of partial viterbi path
  int n = data.nrow(); int dim = data.ncol() - 1;
  double lambda0 = theta[0], lambda1 = theta[1], lambda2 = theta[2];
  double p = theta[4];
  double ps0 = 1. / lambda0 / (1. / lambda0 + p / lambda1 + (1 - p) / lambda2);
  double ps1 = p / lambda1 / (1. / lambda0 + p / lambda1 + (1 - p) / lambda2);
  double ps2 = (1 - p) / lambda2 / (1. / lambda0 + p / lambda1 + (1 - p) / lambda2);
  NumericVector tt = data.column(0);
  NumericMatrix x = data(Range(0, n - 1), Range(1, dim));

  // bf_result matrix: frist three col for forward
  //                   last  three col for backward
  NumericMatrix bf_result(n + 1, 6);
  bf_result(0, 0) = ps0; bf_result(0, 1) = ps1; bf_result(0, 2) = ps2;
  bf_result(n, 3) =   1; bf_result(n, 4) =   1; bf_result(n, 5) =   1;
  NumericVector dx(n);
  
  // result matrix: three cols stand for Viterbi prob of state 0,1,2
  //                at current time points. For numerical reason,
  //                the log-prob is returned.
  NumericMatrix result(pathlength, 3);
  NumericVector cartV(3);
  NumericVector cartW(3);
  
  

  // calculate all h functions
  NumericVector
    hresult00 = ths_h00(x, tt, theta, integrControl),
    hresult01 = ths_h01(x, tt, theta, integrControl),
    hresult02 = ths_h02(x, tt, theta, integrControl),
    hresult10 = ths_h10(x, tt, theta, integrControl),
    hresult11 = ths_h11(x, tt, theta, integrControl),
    hresult12 = ths_h12(x, tt, theta, integrControl),
    hresult20 = ths_h20(x, tt, theta, integrControl),
    hresult21 = ths_h21(x, tt, theta, integrControl),
    hresult22 = ths_h22(x, tt, theta, integrControl);
  
  for (int i = 0; i < n; i++) {
    NumericVector crow = x.row(i);
    if (is_true(all(crow == 0.))) {
      hresult00[i] = 0.;
      hresult01[i] = 0.;
      hresult02[i] = 0.;
      hresult10[i] = 0.;
      hresult11[i] = exp(-lambda1 * tt[i]);
      hresult12[i] = 0.;
      hresult20[i] = 0.;
      hresult21[i] = 0.;
      hresult22[i] = exp(-lambda2 * tt[i]);
    }
  }

  // forward algorithm
  for (int i = 0; i < n; i++) {
    double sumf0 = bf_result(i, 0) * hresult00[i] + bf_result(i, 1) * hresult10[i] + bf_result(i, 2) * hresult20[i];
    double sumf1 = bf_result(i, 0) * hresult01[i] + bf_result(i, 1) * hresult11[i] + bf_result(i, 2) * hresult21[i];
    double sumf2 = bf_result(i, 0) * hresult02[i] + bf_result(i, 1) * hresult12[i] + bf_result(i, 2) * hresult22[i];
    dx[i] = sumf0 + sumf1 + sumf2;
    bf_result(i + 1, 0) = sumf0 / dx[i];
    bf_result(i + 1, 1) = sumf1 / dx[i];
    bf_result(i + 1, 2) = sumf2 / dx[i];
  }

  // backward algorithm
  for (int i = 0; i < n; i++) {
    double sumb0 = bf_result(n-i, 3) * hresult00[n-i-1] + bf_result(n-i, 4) * hresult01[n-i-1] + bf_result(n-i, 5) * hresult02[n-i-1];
    double sumb1 = bf_result(n-i, 3) * hresult10[n-i-1] + bf_result(n-i, 4) * hresult11[n-i-1] + bf_result(n-i, 5) * hresult12[n-i-1];
    double sumb2 = bf_result(n-i, 3) * hresult20[n-i-1] + bf_result(n-i, 4) * hresult21[n-i-1] + bf_result(n-i, 5) * hresult22[n-i-1];
    bf_result(n-i-1, 3) = sumb0 / dx[n-i-1];
    bf_result(n-i-1, 4) = sumb1 / dx[n-i-1];
    bf_result(n-i-1, 5) = sumb2 / dx[n-i-1];
  }

  // prepare for viterbi path
  result(0, 0) = log(bf_result(startpoint, 0));
  result(0, 1) = log(bf_result(startpoint, 1));
  result(0, 2) = log(bf_result(startpoint, 2));
  cartV = result.row(0);
  int ite_stop = startpoint + pathlength - 2;

  // viterbi algorithm
  for (int i = startpoint; i < ite_stop; i++) {
    cartW[0] = cartV[0] + log(hresult00[i]);
    cartW[1] = cartV[1] + log(hresult10[i]);
    cartW[2] = cartV[2] + log(hresult20[i]);
    result(i - startpoint + 1, 0) = max(cartW);

    cartW[0] = cartV[0] + log(hresult01[i]);
    cartW[1] = cartV[1] + log(hresult11[i]);
    cartW[2] = cartV[2] + log(hresult21[i]);
    result(i - startpoint + 1, 1) = max(cartW);
//.........这里部分代码省略.........