Блог

Я пытаюсь решить, имеет ли смысл внедрять оператор R %*% в RCpp, если мой набор данных огромен. НО у меня действительно возникли проблемы с получением реализации RCpp.

Вот мой пример R-кода

 # remove everything in the global environment
rm(list = ls())                

n_states = 4                   # number of states
v_n <- c("H", "S1", "S2", "D") # the 4 states of the model:
n_t = 100                      # number of transitions

# create transition matrix with random numbers. This transition matrix is constant.
m_P = matrix(runif(n_states*n_t), # insert n_states * n_t random numbers (can change this later)
         nrow = n_states, 
         ncol = n_states,
         dimnames = list(v_n, v_n))

# create markov trace, what proportion of population in each state at each period (won't make sense due to random numbers but that is fine)
m_TR <- matrix(NA, 
               nrow = n_t   1 , 
               ncol = n_states, 
               dimnames = list(0:n_t, v_n))     # create Markov trace (n_t   1 because R doesn't understand  Cycle 0)

# initialize Markov trace 
m_TR[1, ] <- c(1, 0, 0, 0)                      


# run the loop

microbenchmark::microbenchmark(   # function from microbenchmark library used to calculate how long this takes to run

for (t in 1:n_t){                              # throughout the number of cycles
  m_TR[t   1, ] <- m_TR[t, ] %*% m_P           # estimate the Markov trace for cycle the next cycle (t   1)
}


) # end of micro-benchmark function

print(m_TR)   # print the result.
  

И вот замена оператора %*%: (Который, похоже, вообще не работает правильно, хотя я не могу понять, почему.

     library(Rcpp)


cppFunction(
    'void estimate_markov(int n_t, NumericMatrix m_P, NumericMatrix m_TR )
     {
            // We want to reproduce this
            //   matrix_A[X 1,]  <- matrix_A[X,] %*% matrix_B
            // The %*% operation behaves as follows for a vector_V %*% matrix_M

            // Here the Matrix M is populated with letters so that you can
            // more easily see how the operation is performed
            // So, a multiplication like this:
            //
            //  V          M
            // {1}  %*%  {A D}
            // {2}       {B E}
            // {3}       {C F}
            //
            // Results in a vector:
            //   V_result
            // {1*A   1*D}
            // {2*B   2*E}
            // {3*C   3*F}
            //
            // Now use values instead of letters for M (ex: A=1, B=2, .., F=6)
            //  V_result
            // {1*1   1*4}    {1   4}     {5}
            // {2*2   2*5} => {4   10} => {14}
            // {3*3   3*6}    {9   18}    {27}
            //
            // Note that the RHS matrix may contain any number of columns,
            // but *MUST* must contain the same number of rows as LHS vector

        // Get dimensions of matricies , and sanity check
        // number of elements in a vector from the LHS matrix must equal == number of row on the RHS
        if( m_TR.cols() != m_P.rows())
            throw std::range_error("Matrix mismatch, m_P.rows != m_TR.rows");

        // we want to know these dimensions, and there is no reason to call these functons in a loop
        // store the values once
        int cnt_P_cols = m_P.cols();
        int cnt_TR_cols = m_TR.cols();

        //
        for(int Index = 1; Index <= n_t;   Index)
        {
            // iterate over the columns in m_TR
            for(int col_iter = 0; col_iter < cnt_TR_cols;   col_iter)
            {
                // an accumulator for the vector multiplication
                double sum = 0;

                // The new value comes from the previous row (Index-1)
                double orig_TR = m_TR(col_iter, Index-1);

                // iterate over the columns in m_P corresponding to this Index
                for(int p_iter = 0; p_iter < cnt_P_cols;   p_iter)
                {
                    // accumulate the value of this TR scalar * the m_P vector
                    sum  = orig_TR * m_P(p_iter, Index);
                }
                m_TR(col_iter, Index) = sum;
            }
          }
        }'
    )
  

Может кто-нибудь указать мне, где моя логика идет не так.