static Rcpp::NumericMatrix x_tilde_3(Rcpp::NumericMatrix X, 
				     Rcpp::IntegerVector nk,
				     Rcpp::IntegerMatrix groups, 
				     Rcpp::NumericMatrix alpha_new,
				     Rcpp::NumericVector d_new)
{
  int K = nk.size();
  int n_tot = X.nrow();
  int p = X.ncol();
  int L = groups.ncol();
  Rcpp::NumericMatrix result(n_tot, L * K);
  
  int idx = 0;
  for (int k = 0; k < K; k++) {
    int n = nk[k];
    for (int l = 0; l < L; l++) {	
      for (int i = 0; i < n; i++) {
	double sum = 0.0;
	for (int j = 0; j < p; j++) {
	  if (elem(groups, j, l)) {
	    sum += elem(X, idx + i, j) * elem(alpha_new, j, k);
	  }
	}
	elem(result, idx + i, L * k + l) = d_new[l] * sum;
      }
    }
    idx += n;
  }	
  return result;
}

void
run_mle(int argc, char* argv[])
{

    // initialize R
    RInside R(argc, argv); 

    // load BradleyTerry library
    // load data to R object
    // Run the MLE
    // Do something with results

    std::string str = 
        "cat('Requireing libraray\n');library('BradleyTerry2'); "
        "cat('Loading data from file\n'); data <- read.table('data',sep=',') ; "
        "cat('Running BTm()\n');fighterModel <- BTm(cbind(win1, win2), fighter1, fighter2, ~ fighter, id='fighter', data=data) ; "
        "BTabilities(fighterModel)"; // returns a matrix of two colums: fighter; ability
    
    Rcpp::NumericMatrix m = R.parseEval(str);   // eval string, return value then as signed to num. vec              

    for (int i=0; i< m.nrow(); i++) { 
        cout << "Figher " << i << " has skill " << m(i,0) << endl;
    }
    cout << endl;

}

//' Marginal correlation matrix
//' 
//' Various workhorse functions to compute the marginal (or unconditional) 
//' correlations (and cross-correlation) estimates efficiently. 
//' They are (almost) 
//' equivalent implementations of \code{\link[stats]{cor}} in Rcpp, 
//' RcppArmadillo, and RcppEigen.
//' 
//' @rdname corFamily
//' @aliases corFamily
//'   corRcpp xcorRcpp corArma xcorArma corEigen xcorEigen
//' @param X A numeric matrix.
//' @param Y A numeric matrix of compatible dimension with the \code{X}, i.e. 
//'   \code{nrow(X)} equals \code{nrow(Y)}.
//' @return
//'   The \code{corXX} family returns a numeric correlation matrix of size 
//'   \code{ncol(X)} times \code{ncol(X)}.
//'   
//'   The \code{xcorXX} family returns a numeric cross-correlation matrix 
//'   of size \code{ncol(X)} times \code{ncol(Y)}.
//' @details
//'   Functions almost like \code{\link{cor}}.
//'   For the \code{xcorXX} functions, the \code{i}'th and \code{j}'th 
//'   entry of the output matrix is the correlation between \code{X[i, ]} and 
//'   \code{X[j, ]}.
//'   Likewise, for the \code{xcorXX} functions, the \code{i}'th and
//'   \code{j}'th entry of the output is the correlation between \code{X[i, ]} 
//'   and \code{Y[j, ]}.
//' @note 
//'   \code{NA}s in \code{X} or \code{Y} will yield \code{NA}s in the correlation matrix.
//'   This also includes the diagonal unlike the behavior of 
//'   \code{\link[stats]{cor}}.
//' @author Anders Ellern Bilgrau <anders.ellern.bilgrau (at) gmail.com>
//' @export
// [[Rcpp::export]]
Rcpp::NumericMatrix corRcpp(Rcpp::NumericMatrix & X) {
  
  const int m = X.ncol();
  const int n = X.nrow();
  
  // Centering the matrix
  X = centerNumericMatrix(X);
  
  Rcpp::NumericMatrix cor(m, m);
  
  // Degenerate case
  if (n == 0) {
    std::fill(cor.begin(), cor.end(), Rcpp::NumericVector::get_na());
    return cor; 
  }
  
  // Compute 1 over the sample standard deviation
  Rcpp::NumericVector inv_sqrt_ss(m);
  for (int i = 0; i < m; ++i) {
    inv_sqrt_ss(i) = 1/sqrt(Rcpp::sum(X(Rcpp::_, i)*X(Rcpp::_, i)));
  }
  
  // Computing the correlation matrix
  for (int i = 0; i < m; ++i) {
    for (int j = 0; j <= i; ++j) {
      cor(i, j) = Rcpp::sum(X(Rcpp::_,i)*X(Rcpp::_,j)) *
        inv_sqrt_ss(i) * inv_sqrt_ss(j);
      cor(j, i) = cor(i, j);
    }
  }

  return cor;
}

//' Check whether there are any non-finite values in a matrix
//'
//' The C++ functions will not work with NA values, and the calculation of the
//' summary profile will take a long time to run before crashing.
//'
//' @param matPtr matrix to check.
//' 
//' @return
//'  Throws an error if any \code{NA}, \code{NaN}, \code{Inf}, or \code{-Inf}
//'  values are found, otherwise returns silently.
//' 
// [[Rcpp::export]]
void CheckFinite(Rcpp::NumericMatrix matPtr) {
  arma::mat mat = arma::mat(matPtr.begin(), matPtr.nrow(), matPtr.ncol(), false, true);
  arma::uvec nonFiniteIdx = arma::find_nonfinite(mat);
  if (nonFiniteIdx.n_elem > 0) {
    throw Rcpp::exception("matrices cannot have non-finite or missing values");
  } 
}

void eigs_sym_shift_c(
    mat_op op, int n, int k, double sigma,
    const spectra_opts *opts, void *data,
    int *nconv, int *niter, int *nops,
    double *evals, double *evecs, int *info
)
{
    BEGIN_RCPP

    CRealShift cmat_op(op, n, data);
    Rcpp::List res;
    try {
        res = run_eigs_shift_sym((RealShift*) &cmat_op, n, k, opts->ncv, opts->rule,
                                 sigma, opts->maxitr, opts->tol, opts->retvec != 0);
        *info = 0;
    } catch(...) {
        *info = 1;  // indicates error
    }

    *nconv = Rcpp::as<int>(res["nconv"]);
    *niter = Rcpp::as<int>(res["niter"]);
    *nops  = Rcpp::as<int>(res["nops"]);
    Rcpp::NumericVector val = res["values"];
    std::copy(val.begin(), val.end(), evals);
    if(opts->retvec != 0)
    {
        Rcpp::NumericMatrix vec = res["vectors"];
        std::copy(vec.begin(), vec.end(), evecs);
    }

    VOID_END_RCPP
}

    /** 
     * Update V, Vtr and fac
     *
     * Note: May want to update fac in a separate operation.  For the
     * fixed-effects modules this will update the factor twice because
     * it is separately updated in updateRzxRx.
     * 
     * @param Xwt square root of the weights for the model matrices
     * @param wtres weighted residuals
     */
    void sPredModule::reweight(Rcpp::NumericMatrix   const&   Xwt,
			       Rcpp::NumericVector   const& wtres) throw(std::runtime_error) {
	if (d_coef.size() == 0) return;
	double one = 1., zero = 0.;
	int Wnc = Xwt.ncol();//, Wnr = Xwt.nrow(),
//	    Xnc = d_X.ncol, Xnr = d_X.nrow;
	if ((Xwt.rows() * Xwt.cols()) != (int)d_X.nrow)
	    throw std::runtime_error("dimension mismatch");
	    // Rf_error("%s: dimension mismatch %s(%d,%d), %s(%d,%d)",
	    // 	     "deFeMod::reweight", "X", Xnr, Xnc,
	    // 	     "Xwt", Wnr, Wnc);
	if (Wnc == 1) {
	    if (d_V) M_cholmod_free_sparse(&d_V, &c);
	    d_V = M_cholmod_copy_sparse(&d_X, &c);
	    chmDn csqrtX(Xwt);
	    M_cholmod_scale(&csqrtX, CHOLMOD_ROW, d_V, &c);
	} else throw runtime_error("sPredModule::reweight: multiple columns in Xwt");
// FIXME write this combination using the triplet representation
	
	chmDn cVtr(d_Vtr);
	const chmDn cwtres(wtres);
	M_cholmod_sdmult(d_V, 'T', &one, &zero, &cwtres, &cVtr, &c);

	CHM_SP Vt = M_cholmod_transpose(d_V, 1/*values*/, &c);
	d_fac.update(*Vt);
	M_cholmod_free_sparse(&Vt, &c);
    }

///********************************************************************
///**  cdm_rcpp_irt_classify_individuals
// [[Rcpp::export]]
Rcpp::List cdm_rcpp_irt_classify_individuals( Rcpp::NumericMatrix like )
{
    int N = like.nrow();
    int TP = like.ncol();
    Rcpp::IntegerVector class_index(N);
    Rcpp::NumericVector class_maxval(N);
    double val=0;
    int ind=0;
    for (int nn=0; nn<N; nn++){
        val=0;
        ind=0;
        for (int tt=0; tt<TP; tt++){
            if ( like(nn,tt) > val ){
                val = like(nn,tt);
                ind = tt;
            }
        }
        class_index[nn] = ind + 1;
        class_maxval[nn] = val;
    }
    //---- OUTPUT:
    return Rcpp::List::create(
                Rcpp::Named("class_index") = class_index,
                Rcpp::Named("class_maxval") = class_maxval
        );
}

static double bic_linear(Rcpp::NumericMatrix X, 
			 Rcpp::NumericVector y, 
			 Rcpp::NumericMatrix beta_new, 
			 double eps, 
			 Rcpp::IntegerVector nk)
{
  int n_tot = X.nrow();
  int p = X.ncol();
  int K = nk.size();
 

  /*calculate SSe*/
  double SSe = 0.0;  
  int idx = 0;
  
  for (int k = 0; k < K; k++) {
    int n = nk[k];
    for (int i = 0; i < n; i++) {
      double Xrow_betacol = 0.0;
      for (int j = 0; j < p; j++) {
	Xrow_betacol += elem(X, idx+i, j) * elem(beta_new, j, k);
      }
      SSe += pow(y[idx+i] - Xrow_betacol, 2);
    }
    idx += n;
  }
  
  double ll = -n_tot / 2.0 * (log(SSe) - log(n_tot) + log(2.0 * M_PI) + 1);
  double bic = -2 * ll + df(beta_new, eps) * log(n_tot);

  return bic;
}

void ScoreGaussL0PenScatter::setData(Rcpp::List& data)
{
	std::vector<int>::iterator vi;
	//uint i;

	// Cast preprocessed data from R list
	dout.level(2) << "Casting preprocessed data...\n";
	_dataCount = Rcpp::as<std::vector<int> >(data["data.count"]);
	dout.level(3) << "# samples per vertex: " << _dataCount << "\n";
	_totalDataCount = Rcpp::as<uint>(data["total.data.count"]);
	dout.level(3) << "Total # samples: " << _totalDataCount << "\n";
	Rcpp::List scatter = data["scatter"];
	Rcpp::NumericMatrix scatterMat;
	_disjointScatterMatrices.resize(scatter.size());
	dout.level(3) << "# disjoint scatter matrices: " << scatter.size() << "\n";
	for (R_len_t i = 0; i < scatter.size(); ++i) {
		scatterMat = Rcpp::NumericMatrix((SEXP)(scatter[i]));
		_disjointScatterMatrices[i] = arma::mat(scatterMat.begin(), scatterMat.nrow(), scatterMat.ncol(), false);
	}

	// Cast index of scatter matrices, adjust R indexing convention to C++
	std::vector<int> scatterIndex = Rcpp::as<std::vector<int> >(data["scatter.index"]);
	for (std::size_t i = 0; i < scatterIndex.size(); ++i)
		_scatterMatrices[i] = &(_disjointScatterMatrices[scatterIndex[i] - 1]);

	// Cast lambda: penalty constant
	_lambda = Rcpp::as<double>(data["lambda"]);
	dout.level(3) << "Penalty parameter lambda: " << _lambda << "\n";

	// Check whether an intercept should be calculated
	_allowIntercept = Rcpp::as<bool>(data["intercept"]);
	dout.level(3) << "Include intercept: " << _allowIntercept << "\n";
}

static double bic_logistic(Rcpp::NumericMatrix X, 
			   Rcpp::NumericVector y, 
			   Rcpp::NumericMatrix beta_new, 
			   double eps, 
			   Rcpp::IntegerVector nk)
{
  int n_tot = X.nrow();
  int p = X.ncol();
  int K = nk.size();
    
  int idx = 0;
  double ll = 0.0;
  for (int k = 0; k < K; k++) {
    int n = nk[k];
    for (int i = 0; i < n; i++) {
      double lp = 0.0;
      for (int j = 0; j < p; j++) {
	lp += elem(X, idx+i, j) * elem(beta_new, j, k);
      }
      ll += y[idx+i] * lp - log(1.0 + exp(lp));
    }
    idx += n;
  }
  
  double bic = -2.0 * ll + df(beta_new, eps) * log(n_tot);
  return bic;
}

///********************************************************************
///** scale2_NA_C
// [[Rcpp::export]]
Rcpp::NumericMatrix scale2_NA_C( Rcpp::NumericMatrix x )
{
    int n = x.nrow();
    int p = x.ncol();
    Rcpp::NumericMatrix y(n,p);
    double mvv=0;
    double sdvv=0;
    double nvv=0;
    double eps_add = 1e-10;
    for (int vv=0;vv<p;vv++){
        mvv=0;
        sdvv=0;
        nvv=0;
        for (int ii=0;ii<n;ii++){
            if (! R_IsNA(x(ii,vv)) ) {
                mvv += x(ii,vv);
                sdvv += std::pow( x(ii,vv), 2.0 );
                nvv ++;
            }
        }
        mvv = mvv / nvv;
        sdvv = std::sqrt( ( sdvv - nvv * mvv*mvv )/(nvv - 1.0 ) );
        // define standardization
        y(_,vv) = ( x(_,vv) - mvv ) / ( sdvv + eps_add );
    }
    //--- output
    return y;
}

//[[Rcpp::export]]
Rcpp::NumericMatrix ZERO(Rcpp::NumericMatrix x) {
  int i=0, j=0;
  for(i=0; i < x.ncol(); i++) {
    for(j=0; j < x.nrow(); j++) {
      x(i,j) = 0;
    }
  }
  return(x);
}

//[[Rcpp::export]]
void decorr(Rcpp::NumericMatrix x) {
  unsigned int i = 1, j=1, n=x.nrow();
  if(n != x.ncol()) Rcpp::stop("matrix is not square");
  for(i=0; i < n; i++) {
    for(j=0; j < n; j++) {
      if(j!=i) x(i,j) = x(i,j)*sqrt(x(i,i)*x(j,j));
    }
  }
}

// kernel Dist function on a Grid
// [[Rcpp::export]]
Rcpp::NumericVector
KdeDist(const Rcpp::NumericMatrix & X
      , const Rcpp::NumericMatrix & Grid
      , const double                h
      , const Rcpp::NumericVector & weight
      , const bool printProgress
	) {
	const unsigned sampleNum = X.nrow();
	const unsigned dimension = Grid.ncol();
	const unsigned gridNum = Grid.nrow();
	// first = sum K_h(X_i, X_j), second = K_h(x, x), third = sum K_h(x, X_i)
	std::vector< double > firstValue;
	const double second = 1.0;
	std::vector< double > thirdValue;
	double firstmean;
	Rcpp::NumericVector kdeDistValue(gridNum);
	int counter = 0, percentageFloor = 0;
	int totalCount = sampleNum + gridNum;

	if (printProgress) {
		printProgressFrame(Rprintf);
	}

	firstValue = computeKernel< std::vector< double > >(
			X, X, h, weight, printProgress, Rprintf, counter, totalCount,
			percentageFloor);

	if (dimension <= 1) {
		thirdValue = computeKernel< std::vector< double > >(
				X, Grid, h, weight, printProgress, Rprintf, counter, totalCount,
				percentageFloor);
	}
	else {
		thirdValue = computeGaussOuter< std::vector< double > >(
				X, Grid, h, weight, printProgress, Rprintf, counter, totalCount,
				percentageFloor);
	}

	if (weight.size() == 1) {
		firstmean = std::accumulate(firstValue.begin(), firstValue.end(), 0.0) / sampleNum;
	}
	else {
		firstmean = std::inner_product(
				firstValue.begin(), firstValue.end(), weight.begin(), 0.0) / 
				std::accumulate(weight.begin(), weight.end(), 0.0);
	}

	for (unsigned gridIdx = 0; gridIdx < gridNum; ++gridIdx) {
		kdeDistValue[gridIdx] = std::sqrt(firstmean + second - 2 * thirdValue[gridIdx]);
	}

	if (printProgress) {
		Rprintf("\n");
	}

	return kdeDistValue;
}

static void print(Rcpp::NumericMatrix A)
{
  for (int i = 0; i < A.nrow(); i++) {
    for (int j = 0; j < A.ncol(); j++) {
      printf("%9f ", elem(A, i, j));
    }
    putchar('\n');
  }
}

static int df(Rcpp::NumericMatrix beta_new, double eps)
{
  int result = 0;
  int n = beta_new.nrow() * beta_new.ncol();
  for(int i=0; i < n; i++){
    result += nz(beta_new[i],eps);
  }
  return result;
}

// [[Rcpp::export]]
Rcpp::NumericMatrix testColPost(Rcpp::NumericMatrix post, Rcpp::List m2u, int nthreads){
    Rcpp::IntegerVector values = Rcpp::as<Rcpp::IntegerVector>(m2u["values"]);
    Rcpp::IntegerVector map = Rcpp::as<Rcpp::IntegerVector>(m2u["map"]);
    if (post.ncol() != map.length()) Rcpp::stop("posteriors doesn't match with m2u");
    
    Rcpp::NumericMatrix smallerPost(post.nrow(), values.length());
    Vec<double> foo; NMPreproc preproc(asVec(values), asVec(map), foo);
    collapsePosteriors_core(asMat(smallerPost), asMat(post), preproc);
    return smallerPost;
}

// gradient-descent ------------------------------------------------------------
//' @title Regular gradient descent for latent factor model with covariates
//' @param X The ratings matrix. Unobserved entries must be marked NA. Users
//' must be along rows, and tracks must be along columns.
//' @param Z The covariates associated with each pair. This must be shaped as an
//' array with users along rows, tracks along columns, and features along
//' slices.
//' @param k_factors The number of latent factors for the problem.
//' @param lambdas The regularization parameters for P, Q, and beta, respectively.
//' @param n_iter The number of gradient descent iterations to run.
//' @param batch_samples The proportion of the observed entries to sample in the
//' gradient for each factor in the gradient descent.
//' @param batch_factors The proportion of latent factors to update at each
//' iteration.
//' @param gamma The step-size in the gradient descent.
//' @param verbose Print the objective at each iteration of the descent?
//' @return A list with the following elements, \cr
//'   $P The learned user latent factors. \cr
//'   $Q The learned track latent factors. \cr
//'   $beta The learned regression coefficients.
//' @export
// [[Rcpp::export]]
Rcpp::List svd_cov(Rcpp::NumericMatrix X, Rcpp::NumericVector Z_vec,
		   int k_factors, Rcpp::NumericVector lambdas, int n_iter,
		   double batch_samples, double batch_factors,
		   double gamma_pq, double gamma_beta, bool verbose) {

  // convert to arma
  Rcpp::IntegerVector Z_dim = Z_vec.attr("dim");
  arma::cube Z(Z_vec.begin(), Z_dim[0], Z_dim[1], Z_dim[2], false);

  // initialize results
  int m = X.nrow();
  int n = X.ncol();
  int l = Z_dim[2];

  arma::mat P = 1.0 / sqrt(m) * arma::randn(m, k_factors);
  arma::mat Q = 1.0 / sqrt(n) * arma::randn(n, k_factors);
  arma::vec beta = 1.0 / sqrt(l) * arma::randn(l);
  arma::vec objs = arma::zeros(n_iter);

  // perform gradient descent steps
  for(int cur_iter = 0; cur_iter < n_iter; cur_iter++) {
    // update user factors
    arma::uvec cur_factors = get_batch_ix(m, batch_factors);
    for(int i = 0; i < cur_factors.n_elem; i++) {
      int cur_ix = cur_factors(i);
      P.row(cur_ix) = update_factor(X.row(cur_ix),
				    Z(arma::span(cur_ix), arma::span(), arma::span()),
				    arma::conv_to<arma::vec>::from(P.row(cur_ix)), Q, beta,
				    batch_samples, lambdas[0], gamma_pq).t();
    }

    // update track factors
    cur_factors = get_batch_ix(n, batch_factors);
    for(int j = 0; j < cur_factors.n_elem; j++) {
      int cur_ix = cur_factors(j);
      Q.row(cur_ix) = update_factor(X.column(cur_ix),
				    Z(arma::span(), arma::span(cur_ix), arma::span()),
				    arma::conv_to<arma::vec>::from(Q.row(cur_ix)), P, beta,
				    batch_samples, lambdas[1], gamma_pq).t();
    }

    // update regression coefficients
    beta = update_beta(Rcpp::as<arma::mat>(X), Z, P, Q, beta, lambdas[2], gamma_beta);
    objs[cur_iter] = objective_fun(Rcpp::as<arma::mat>(X), Z, P, Q, beta, lambdas);
    if(verbose) {
      printf("iteration %d \n", cur_iter);
      printf("Objective: %g \n", objs[cur_iter]);
    }
  }

  return Rcpp::List::create(Rcpp::Named("P") = P,
			    Rcpp::Named("Q") = Q,
			    Rcpp::Named("beta") = beta,
			    Rcpp::Named("objs") = objs);
}

////// Ordinary DistanceMapper //////////////
R_DimSqueezer::R_DimSqueezer(Rcpp::NumericMatrix r_positions)
{
  mapper = 0;
  multithreaded = true;
  node_no = r_positions.nrow();
  dimension_no = r_positions.ncol();
  positions = matrix_to_array(r_positions);
  distances = inter_node_distances(positions, node_no, dimension_no);

  if(positions)
    mapper = new DistanceMapper(node_no, dimension_no, positions, distances);
}

void sampleRho(Rcpp::NumericVector& rho, Rcpp::NumericMatrix& A, arma::mat A_restrict, double rhoa, double rhob){
    int k = A.ncol();
    int p = A.nrow();
	arma::rowvec A_maxnnz_col = arma::sum(A_restrict, 0);
    for (int i=0; i<k; i++) {
        NumericMatrix::Column col = A.column(i) ;
        NumericVector As = ifelse(abs(col)<1e-10, 0.0, 1.0);
        double nnz = std::accumulate(As.begin(), As.end(), 0.0);
        double maxnnz = A_maxnnz_col(i);
        rho(i) = Rf_rbeta(rhoa + nnz, rhob + maxnnz-nnz);
    }
}