Hi everyone,

I've just heard about Eigen recently and it looks really promising,therefore I decided to used it for our program (a bioinformatics tool) by transforming the performance-critical codes (our program could need a week to compute some datasets) into Eigen's notations.

I've spent sometimes thinking about an elegant solution but still cannot come up with any reasonable idea, probably due to lacking of experience with Eigen and linear algebra in general. Here is a snippet of the code inside a function (which is called 10000s times):

Code: Select all

   double tree_lh = node_branch->lh_scale_factor + dad_branch->lh_scale_factor;
   int ncat = site_rate->getNRate();
   double p_invar = site_rate->getPInvar();
   double p_var_cat = (1.0 - p_invar) / (double)ncat;
   int nstates = aln->num_states;
   int block = ncat * nstates;
   int trans_size = nstates * nstates;
   int ptn, cat, state1, state2;

   double trans_mat[ncat*trans_size];
   double state_freq[nstates];
   model->getStateFrequency(state_freq);

   for (cat = 0; cat < ncat; cat++) {
      //trans_mat[cat] = model->newTransMatrix();
      double *trans_cat = trans_mat + (cat*trans_size);
      model_factory->computeTransMatrix(dad_branch->length * site_rate->getRate(cat), trans_cat);
      for (state1 = 0; state1 < nstates; state1++) {
         double *trans_mat_state = trans_cat + (state1*nstates);
         for (state2 = 0; state2 < nstates; state2++)
            trans_mat_state[state2] *= state_freq[state1];
      }
   }

   for (ptn = 0; ptn < aln->size(); ptn++) {
      double lh_ptn = 0.0; // likelihood of the pattern
      int dad_state = 1000; // just something big enough
      if (dad->isLeaf())
         dad_state = (*aln)[ptn][dad->id];
      for (cat = 0; cat < ncat; cat++) {
         double *partial_lh_site = node_branch->partial_lh + (ptn*block + cat*nstates);
         double *partial_lh_child = dad_branch->partial_lh + (ptn*block + cat*nstates);
         if (dad_state  < nstates) { // single state
            // external node
            double *trans_state = trans_mat + (cat*trans_size + dad_state*nstates);
            for (state2 = 0; state2 < nstates; state2++)
               lh_ptn += partial_lh_child[state2] * trans_state[state2];
         } else {
            // internal node, or external node but ambiguous character
            for (state1 = 0; state1 < nstates; state1++) {
               double lh_state = 0.0;  // likelihood of state1
               double *trans_state = trans_mat + (cat*trans_size + state1*nstates);
   
               for (state2 = 0; state2 < nstates; state2++)
                  lh_state += partial_lh_child[state2] * trans_state[state2];
                                       
               lh_ptn += lh_state * partial_lh_site[state1];
            }
         }
      }
      lh_ptn *= p_var_cat;
      if ((*aln)[ptn].is_const && (*aln)[ptn][0] < nstates) {
         lh_ptn += p_invar * state_freq[(*aln)[ptn][0]];
      }
      assert(lh_ptn > 0);
      tree_lh += log(lh_ptn) * (*aln)[ptn].frequency;
   }
   return tree_lh;


The code is relatively straightforward : I have some big arrays, which were already initialized beforehand (in the first for loop). The performance-critical code is the second loop for (ptn = 0; ptn < aln->size(); ptn++), the loop size could be 10000s, depend on the dataset. Inside this loop there is a second innerloop for (cat = 0; cat < ncat; cat++), ncat is usually 4.

Until now I am only able to convert code inside this innerloop into Eigen's matrix operations,which is kind of obvious(using 4x4 or 20x20 matrices since nstates is either 4 or 20)

I want to also convert the 2 outerloops ( for (cat..) and for(ptn...) ) into Eigen, but still have no clue !

I would be very nice if someone can help me out

Thanks

The code probably looks scary at first, but only the loop for (ptn = 0; ptn < aln->size(); ptn++) is relevant, and what I want to compute in each loop is the value store in the variable lh_ptn and after this loop ends, the values of lh_ptn are summed up in tree_lh.

If you have time, please spend 5 minutes looking through that loop and tell me what options in Eigen do I have ?

Thanks

This is a computation for a maximum likelihood method, right?

I'd say pull the if's out of the loops and you have two loops left, one matrix/matrix (reduced) and one matrix/vector.

Yes, it is exactly the computation for maximum likelihood. You meant the if (dad_state < nstates) statement right ? I already noticed that it is inefficient and pulled it out in the current code, I didnt want to post the new code here because it is longer, so just ignore the if, just consider the else branch, which is:

Code: Select all

            // internal node, or external node but ambiguous character
            for (state1 = 0; state1 < nstates; state1++) {
               double lh_state = 0.0;  // likelihood of state1
               double *trans_state = trans_mat + (cat*trans_size + state1*nstates);
   
               for (state2 = 0; state2 < nstates; state2++)
                  lh_state += partial_lh_child[state2] * trans_state[state2];
                                       
               lh_ptn += lh_state * partial_lh_site[state1];
            }


I had no problem converting the above code into Eigen matrixs. But this code is inside 2 for loops, which simply kind of sum up all the results computed in the code above. This 2 loops have altogether the size of cat*aln->size(), which is pretty big, cat is usually 4, but aln->size() can be 1000s.

What I wanted to ask is, whether it is possible or meaningful if I try to convert these 2 loops into Eigen Matrixs so that they can be vectorize, loop unrolled, etc... or just let it be I hope the compiler is smart enough to vectorize the code ?

Thanks

I'm not sure (you should just try it out), but I feel like it won't buy you much to transform the outer-most loop into matrix operations. As the code is already out of date, I can't tell you specific things, but as I said the code you posted first has a lot of opportunity for classic optimisation.

I think I have managed to come up with a reasonable solution/optimization.

@Andre: What do you mean by classic optimization ?