[Debian-med-packaging] MP patch for PhyML

Mon Mar 14 19:11:38 UTC 2011

Hi Andreas,

No, I haven't had the time to test it yet.
But I will look into it very soon.

Regards,

-Stephane-

On 15/03/11 07:02, Andreas Tille wrote:
> Hi,
> 
> I noticed that there is a new version of PhyML (20110304) available.
> Does this regard the MP feature yet?  (I have not checked the source
> myself - that's why I'm asking.)
> 
> Kind regards
> 
>      Andreas.
> 
> On Tue, Feb 15, 2011 at 11:38:16PM +0100, Diego Darriba López wrote:
>> Dear Andreas,
>>
>>
>> My OpenMP implementation perform the MP just using pragma directives, so if the -fopenmp flag is not set, the compiler will just ignore these lines and it will perform a single processor build. I checked both builds with gcc. We would be in a different situation if I had used OpenMP API calls. In fact, the #include<omp.h> I included in my patch it is not necessary at all (sorry). This way, the MP or sequential build just depends on the use of the -fopenmp flag.
>>
>>
>> Keep in mind that this flag differs between compilares (e.g, -openmp for icc, -mp for openUH, etc.) if you want to bring support for different compilers in the autoconf file.
>>
>>
>> Kind regards,
>>
>>
>> Diego.
>>
>>
>> De: "Andreas Tille" <andreas at an3as.eu>
>> Para: "Diego Darriba" <ddarriba at udc.es>
>> CC: "s guindon" <s.guindon at auckland.ac.nz>, "Stephane Guindon" <guindon at stat.auckland.ac.nz>, "Olivier Gascuel" <gascuel at lirmm.fr>, "Jean-Francois Dufayard" <jeanfrancois.dufayard at gmail.com>, "David Posada" <dposada at uvigo.es>, "Debian Med Packaging Team" <debian-med-packaging at lists.alioth.debian.org>
>> Enviados: Martes, 15 de Febrero 2011 21:25:51
>> Asunto: Re: MP patch for PhyML
>>
>> Dear Diego,
>>
>> I admit I have no experience with MP builds but I wonder whether you
>> should not include the changes into #ifdef USEMP / #else / #endif
>> statements to enable both - single processor and MP builds. Or does the
>> patch work with both and it is just the build option which controls
>> everything?
>>
>> Kind regards
>>
>> Andreas.
>>
>> On Tue, Feb 15, 2011 at 06:58:19PM +0100, Diego Darriba wrote:
>>> Dear Stephan et al,
>>>
>>> I send you the patch for the current phyml trunk. I checked the results
>>> with several input data, and I think I solved correctly some
>>> dependencies in the parallelized loops from my computer science point of
>>> view. In the configure.ac it is just necessary to include the *-fopenmp*
>>> flag for compiling (if using gcc compiler). I would modify it, but I
>>> prefer you to do so. That way you can create a conditional Makefile, or
>>> include the flag as a commentary, or whatever.
>>>
>>> If the compilation flag is not set, the application will compile with no
>>> error, but it will run as before (i.e, sequentially).
>>>
>>> Please, let me know if you have any doubt or any problem.
>>>
>>> Kind Regards,
>>>
>>> Diego.
>>>
>>> El 15/02/11 02:51, Stephane Guindon escribió:
>>>> Dear Andreas,
>>>>
>>>> Thank you very much for contributing to PhyML.
>>>> Yes, it makes perfect sense to try and include this functionality
>>>> into the latest release of the program. Please feel free to take
>>>> a look at the most recent version of the sources on our svn server
>>>> (http://code.google.com/p/phyml/source/browse/#svn%2Ftrunk%2Fsrc)
>>>> and let me know how to amend the configure.ac, *.c and *.h file
>>>> in an appropriate manner.
>>>>
>>>> Regards,
>>>>
>>>> -Stephane-
>>>>
>>>> On 15/02/11 03:26, Andreas Tille wrote:
>>>>> Hi PhyML authors,
>>>>>
>>>>> when discssing with Diego about a program which uses PhyML he told me
>>>>> about his MP patch which would enable making use of hybrid memory when
>>>>> running PhyML on clusters. I was considering to port this patch to the
>>>>> Debian packaged version where I intended to support two versions of
>>>>> the package: one conventional phylip and one phylip-mp.
>>>>>
>>>>> When inspecting the patch I noticed that it was done against an older
>>>>> version of PhyML (20100123). I would consider it a reasonable idea if
>>>>> this patch would be included into the upstream version by just wrapping
>>>>> it into "#ifdef"s properly to enable switching between both versions by
>>>>> simply providing different options to make. If the patch would be taken
>>>>> over into the official PhyML code this would simplify our work as
>>>>> distributor (as well as Diego's work) and it might be of extra value for
>>>>> other users.
>>>>>
>>>>> What do you think about this idea?
>>>>>
>>>>> Kind regards
>>>>>
>>>>> Andreas.
>>>>>
>>>
>>
>>> Index: src/lk.c
>>> ===================================================================
>>> --- src/lk.c (revisi??n: 582)
>>> +++ src/lk.c (copia de trabajo)
>>> @@ -11,8 +11,8 @@
>>> */
>>>
>>> #include "lk.h"
>>> +#include <omp.h>
>>>
>>> -
>>> /* int LIM_SCALE; */
>>> /* phydbl LIM_SCALE_VAL; */
>>> /* phydbl BIG; */
>>> @@ -320,8 +320,9 @@
>>>
>>> phydbl Lk(t_tree *tree)
>>> {
>>> - int br;
>>> + int br, site;
>>> int n_patterns;
>>> + phydbl c_lnL;
>>>
>>> tree->old_lnL = tree->c_lnL;
>>>
>>> @@ -349,10 +350,17 @@
>>>
>>> tree->c_lnL = .0;
>>> tree->sum_min_sum_scale = .0;
>>> - For(tree->curr_site,n_patterns)
>>> +
>>> + c_lnL = .0;
>>> + #pragma omp parallel for reduction(+:c_lnL)
>>> + For(site,n_patterns)
>>> {
>>> - if(tree->data->wght[tree->curr_site] > SMALL) Lk_Core(tree->noeud[0]->b[0],tree);
>>> + if(tree->data->wght[site] > SMALL) {
>>> + Lk_Core(tree->noeud[0]->b[0],tree,site);
>>> + c_lnL += tree->c_lnL_sorted[site];
>>> + }
>>> }
>>> + tree->c_lnL = c_lnL;
>>>
>>> /* Qksort(tree->c_lnL_sorted,NULL,0,n_patterns-1); */
>>>
>>> @@ -372,7 +380,8 @@
>>>
>>> phydbl Lk_At_Given_Edge(t_edge *b_fcus, t_tree *tree)
>>> {
>>> - int n_patterns;
>>> + int n_patterns,site;
>>> + phydbl c_lnL;
>>>
>>> n_patterns = tree->n_pattern;
>>>
>>> @@ -396,10 +405,17 @@
>>>
>>> tree->c_lnL = .0;
>>> tree->sum_min_sum_scale = .0;
>>> - For(tree->curr_site,n_patterns)
>>> +
>>> + c_lnL = .0;
>>> + #pragma omp parallel for reduction(+:c_lnL)
>>> + For(site,n_patterns)
>>> {
>>> - if(tree->data->wght[tree->curr_site] > SMALL) Lk_Core(b_fcus,tree);
>>> + if(tree->data->wght[site] > SMALL) {
>>> + Lk_Core(b_fcus,tree,site);
>>> + c_lnL += tree->c_lnL_sorted[site];
>>> + }
>>> }
>>> + tree->c_lnL = c_lnL;
>>>
>>> /* Qksort(tree->c_lnL_sorted,NULL,0,n_patterns-1); */
>>>
>>> @@ -422,7 +438,7 @@
>>>
>>> #ifndef USE_OLD_LK
>>>
>>> -phydbl Lk_Core(t_edge *b, t_tree *tree)
>>> +phydbl Lk_Core(t_edge *b, t_tree *tree, int site)
>>> {
>>> phydbl log_site_lk;
>>> phydbl site_lk_cat, site_lk;
>>> @@ -431,7 +447,7 @@
>>> phydbl max_sum_scale;
>>> phydbl sum;
>>> int ambiguity_check,state;
>>> - int catg,ns,k,l,site;
>>> + int catg,ns,k,l;
>>> int dim1,dim2,dim3;
>>> int *sum_scale_left_cat,*sum_scale_rght_cat;
>>> phydbl multiplier;
>>> @@ -439,6 +455,7 @@
>>> phydbl tmp;
>>> phydbl logbig;
>>> phydbl inv_site_lk;
>>> + phydbl tree_site_lk_cat[tree->mod->n_catg];
>>>
>>> logbig = LOG((phydbl)BIG);
>>>
>>> @@ -451,7 +468,6 @@
>>> site_lk_cat = .0;
>>> ambiguity_check = -1;
>>> state = -1;
>>> - site = tree->curr_site;
>>> ns = tree->mod->ns;
>>>
>>> if((b->rght->tax) && (!tree->mod->s_opt->greedy))
>>> @@ -537,7 +553,7 @@
>>> }
>>> }
>>> }
>>> - tree->site_lk_cat[catg] = site_lk_cat;
>>> + tree_site_lk_cat[catg] = site_lk_cat;
>>> }
>>>
>>> max_sum_scale = (phydbl)BIG;
>>> @@ -562,7 +578,7 @@
>>> Warn_And_Exit("\n");
>>> }
>>>
>>> - tmp = sum + (logbig - LOG(tree->site_lk_cat[catg]))/(phydbl)LOG2;
>>> + tmp = sum + (logbig - LOG(tree_site_lk_cat[catg]))/(phydbl)LOG2;
>>> if(tmp < max_sum_scale) max_sum_scale = tmp; /* min of the maxs */
>>> }
>>>
>>> @@ -575,7 +591,7 @@
>>> {
>>> exponent = -(sum_scale_left_cat[catg]+sum_scale_rght_cat[catg])+fact_sum_scale;
>>>
>>> - site_lk_cat = tree->site_lk_cat[catg];
>>> + site_lk_cat = tree_site_lk_cat[catg];
>>>
>>> if(exponent > 0)
>>> {
>>> @@ -605,12 +621,12 @@
>>> PhyML_Printf("\n+ site=%4d cat=%4d site_lk_cat=%G sum_scale=%d max=%G fact=%d expo=%d dbl=%G",
>>> tree->curr_site,
>>> catg,
>>> - tree->site_lk_cat[catg],
>>> + tree_site_lk_cat[catg],
>>> sum_scale_left_cat[catg]+sum_scale_rght_cat[catg],
>>> max_sum_scale,
>>> fact_sum_scale,
>>> -(sum_scale_left_cat[catg]+sum_scale_rght_cat[catg])+fact_sum_scale,
>>> - (double)tree->site_lk_cat[catg] * pow(2.,-(sum_scale_left_cat[catg]+sum_scale_rght_cat[catg])+fact_sum_scale));
>>> + (double)tree_site_lk_cat[catg] * pow(2.,-(sum_scale_left_cat[catg]+sum_scale_rght_cat[catg])+fact_sum_scale));
>>>
>>> site_lk_cat = BIG / 10;
>>> }
>>> @@ -623,12 +639,12 @@
>>> PhyML_Printf("\n- site=%4d cat=%4d site_lk_cat=%G sum_scale=%d max=%G fact=%d expo=%d dbl=%G",
>>> tree->curr_site,
>>> catg,
>>> - tree->site_lk_cat[catg],
>>> + tree_site_lk_cat[catg],
>>> sum_scale_left_cat[catg]+sum_scale_rght_cat[catg],
>>> max_sum_scale,
>>> fact_sum_scale,
>>> -(sum_scale_left_cat[catg]+sum_scale_rght_cat[catg])+fact_sum_scale,
>>> - (double)tree->site_lk_cat[catg] * pow(2.,-(sum_scale_left_cat[catg]+sum_scale_rght_cat[catg])+fact_sum_scale));
>>> + (double)tree_site_lk_cat[catg] * pow(2.,-(sum_scale_left_cat[catg]+sum_scale_rght_cat[catg])+fact_sum_scale));
>>>
>>> Exit("\n");
>>> }
>>> @@ -636,13 +652,13 @@
>>> site_lk_cat = .0;
>>> }
>>>
>>> - tree->site_lk_cat[catg] = site_lk_cat;
>>> + tree_site_lk_cat[catg] = site_lk_cat;
>>> }
>>>
>>> site_lk = .0;
>>> For(catg,tree->mod->n_catg)
>>> {
>>> - site_lk += tree->site_lk_cat[catg] * tree->mod->gamma_r_proba[catg];
>>> + site_lk += tree_site_lk_cat[catg] * tree->mod->gamma_r_proba[catg];
>>> }
>>>
>>> /* The substitution model does include invariable sites */
>>> @@ -677,7 +693,7 @@
>>>
>>> log_site_lk = LOG(site_lk) - (phydbl)LOG2 * fact_sum_scale;
>>>
>>> - For(catg,tree->mod->n_catg) tree->log_site_lk_cat[catg][site] = LOG(tree->site_lk_cat[catg]) - (phydbl)LOG2 * fact_sum_scale;
>>> + For(catg,tree->mod->n_catg) tree->log_site_lk_cat[catg][site] = LOG(tree_site_lk_cat[catg]) - (phydbl)LOG2 * fact_sum_scale;
>>>
>>> if(isinf(log_site_lk) || isnan(log_site_lk))
>>> {
>>> @@ -694,7 +710,6 @@
>>> /* Multiply log likelihood by the number of times this site pattern is found in the data */
>>> tree->c_lnL_sorted[site] = tree->data->wght[site]*log_site_lk;
>>>
>>> - tree->c_lnL += tree->data->wght[site]*log_site_lk;
>>> /* tree->sum_min_sum_scale += (int)tree->data->wght[site]*min_sum_scale; */
>>>
>>> return log_site_lk;
>>> @@ -1118,6 +1133,7 @@
>>> Pij2 = d->b[dir2]->Pij_rr;
>>>
>>> /* For every site in the alignment */
>>> + #pragma omp parallel for private(state_v1,state_v2,ambiguity_check_v1,ambiguity_check_v2,catg,i,p1_lk1,p2_lk2)
>>> For(site,n_patterns)
>>> {
>>> state_v1 = state_v2 = -1;
>>> @@ -1399,6 +1415,7 @@
>>> Pij2 = d->b[dir2]->Pij_rr;
>>>
>>> /* For every site in the alignment */
>>> + #pragma omp parallel for private(state_v1,state_v2,ambiguity_check_v1,ambiguity_check_v2,catg,i,p1_lk1,p2_lk2)
>>> For(site,n_patterns)
>>> {
>>> state_v1 = state_v2 = -1;
>>> @@ -2234,13 +2251,13 @@
>>> sides of t_edge *b are up-to-date. Calculate the log-likelihood at one site.
>>> */
>>> #ifdef USE_OLD_LK
>>> -phydbl Lk_Core(t_edge *b, t_tree *tree)
>>> +phydbl Lk_Core(t_edge *b, t_tree *tree, int site)
>>> {
>>> phydbl log_site_lk, site_lk, site_lk_cat;
>>> phydbl scale_left, scale_rght;
>>> phydbl sum;
>>> int ambiguity_check,state;
>>> - int catg,ns,k,l,site;
>>> + int catg,ns,k,l;
>>> int dim1,dim2,dim3;
>>>
>>>
>>> @@ -2253,7 +2270,6 @@
>>> site_lk_cat = .0;
>>> ambiguity_check = -1;
>>> state = -1;
>>> - site = tree->curr_site;
>>> ns = tree->mod->ns;
>>>
>>> tree->old_site_lk[site] = tree->cur_site_lk[site];
>>> @@ -2383,7 +2399,7 @@
>>> tree->cur_site_lk[site] = log_site_lk;
>>> /* Multiply log likelihood by the number of times this site pattern is found oin the data */
>>> tree->c_lnL_sorted[site] = tree->data->wght[site]*log_site_lk;
>>> - tree->c_lnL += tree->c_lnL_sorted[site];
>>> +
>>> return log_site_lk;
>>> }
>>>
>>> @@ -2492,6 +2508,7 @@
>>> Pij2 = d->b[dir2]->Pij_rr;
>>>
>>> /* For every site in the alignment */
>>> + #pragma omp parallel for private(scale_v1,scale_v2,state_v1,state_v2,ambiguity_check_v1,ambiguity_check_v2,catg,i,j,p1_lk1,p2_lk2,max_p_lk)
>>> For(site,n_patterns)
>>> {
>>> /* Current scaling values at that site */
>>> Index: src/lk.h
>>> ===================================================================
>>> --- src/lk.h (revisi??n: 582)
>>> +++ src/lk.h (copia de trabajo)
>>> @@ -51,7 +51,7 @@
>>> phydbl Update_Lk_At_Given_Edge(t_edge *b_fcus, t_tree *tree);
>>> void Update_P_Lk_Greedy(t_tree *tree, t_edge *b_fcus, t_node *n);
>>> void Get_All_Partial_Lk_Scale_Greedy(t_tree *tree, t_edge *b_fcus, t_node *a, t_node *d);
>>> -phydbl Lk_Core(t_edge *b, t_tree *tree);
>>> +phydbl Lk_Core(t_edge *b, t_tree *tree, int site);
>>> phydbl Lk_Triplet(t_node *a, t_node *d, t_tree *tree);
>>> void Print_Lk_Given_Edge_Recurr(t_node *a, t_node *d, t_edge *b, t_tree *tree);
>>> phydbl *Post_Prob_Rates_At_Given_Edge(t_edge *b, phydbl *post_prob, t_tree *tree);
>>
>>
>> --
>> http://fam-tille.de
> 
> --
> http://fam-tille.de

-- 
Stephane Guindon
Department of Statistics
University of Auckland
http://www.stat.auckland.ac.nz