4 * (c) 2003 Gangolf Jobb and Korbinian Strimmer
6 * Functions for nonparametric rate smoothing (NPRS)
7 * (see MJ Sanderson. 1997. MBE 14:1218-1231)
9 * This code may be distributed under the GNU GPL
17 /*============================================================================*/
24 #define MAXPARAM 1-EPS
25 #define STARTPARAM 0.5
29 #define MINLPARAM log(MINPARAM)
30 #define MAXLPARAM log(MAXPARAM)
31 #define STARTLPARAM log(STARTPARAM)
36 #define index aperates_index
38 int tree_lowerNodes[ARRLEN];
39 int tree_upperNodes[ARRLEN];
40 double tree_edgeLengths[ARRLEN];
42 char tree_tipLabels[ARRLEN][LABLEN];
45 int nparams,index[ARRLEN],ancestor[ARRLEN];
47 /*============================================================================*/
53 double *minEdgeLength,
63 for(i=0;i<tree_nedges;i++) {
64 tree_lowerNodes[i]=lowerNodes[i];
65 tree_upperNodes[i]=upperNodes[i];
66 tree_edgeLengths[i]=(edgeLengths[i]<*minEdgeLength)?*minEdgeLength:edgeLengths[i];
71 for(i=0;i<tree_ntips;i++) {
72 strcpy(tree_tipLabels[i],tipLabels[i]);
76 for(i=0;i<tree_nedges;i++) {
77 if(tree_lowerNodes[i]<0&&tree_upperNodes[i]<0) {index[-tree_upperNodes[i]]=nparams++;}
78 if(tree_upperNodes[i]<0) ancestor[-tree_upperNodes[i]]=tree_lowerNodes[i];
84 /*============================================================================*/
86 void getNFreeParams(int *result) { *result=nparams; }
88 /*============================================================================*/
90 void getNEdges(int *result) { *result=tree_nedges; }
92 /*============================================================================*/
94 void getEdgeLengths(double *result) {
97 for(i=0;i<tree_nedges;i++) result[i]=tree_edgeLengths[i];
101 /*============================================================================*/
103 double age(double *params,int node) {
106 if(node>=0) return(0.0);
107 if(node==-1) return(1.0);
110 while(node!=-1) {prod+=params[index[-node]]; node=ancestor[-node];}
116 void getDurations(double *params,double *scale,double *result) {
119 for(i=0;i<tree_nedges;i++) {
120 low=tree_lowerNodes[i]; upp=tree_upperNodes[i];
121 if(low<0&&upp<0) result[i]=*scale*(age(params,low)-age(params,upp));
122 else result[i]=*scale*age(params,low);
127 /*============================================================================*/
129 void getRates(double *params,double *scale,double *result) {
132 for(i=0;i<tree_nedges;i++) {
133 low=tree_lowerNodes[i]; upp=tree_upperNodes[i];
134 if(low<0&&upp<0) result[i]=tree_edgeLengths[i]/(*scale*(age(params,low)-age(params,upp)));
135 else result[i]=tree_edgeLengths[i]/(*scale*age(params,low));
140 /*============================================================================*/
143 /* note on the choice of parameters:
145 * in order to obtain a chronogram we optimize the
146 * relative heights of each internal node, i.e. to
147 * each internal node (minus the root) we assign a number
148 * between 0 and 1 (MINPARAM - MAXPARAM).
150 * To obtain the actual height of a node the relative heights
151 * of the nodes have to multiplied (in fact we work on the log-scale
156 /* internal objective function - parameters are on log scale */
157 void nprsObjectiveFunction(
163 double me,rk,rj,sum,scale=1.0,durations[ARRLEN],rates[ARRLEN];
167 getDurations(params,&scale,durations);
169 for(k=0;k<tree_nedges;k++) rates[k]=tree_edgeLengths[k]/durations[k];
171 me=0.; for(k=0;k<tree_nedges;k++) me+=rates[k]; me/=tree_nedges;
175 for(k=0;k<tree_nedges;k++) { rk=rates[k];
176 for(j=0;j<tree_nedges;j++) { if(j==k) continue; rj=rates[j];
177 if(tree_lowerNodes[j]==-1) sum+=pow(fabs(me-rj),p); else
178 if(tree_upperNodes[k]==tree_lowerNodes[j]) sum+=pow(fabs(rk-rj),p);
187 /* check parameter bounds on log scale */
188 int checkLogParams(double *params)
191 for(i=0; i<nparams; i++)
193 if(params[i] > MAXLPARAM || params[i] < MINLPARAM) return 0; /* out of bounds */
195 return 1; /* within bounds */
201 * public objective function
202 * - parameters are on log scale
203 * - if parameters are out of bounds function returns large value
205 void objFuncLogScale(
212 if( checkLogParams(params) == 0 ) /* out of bounds */
218 nprsObjectiveFunction(params, expo, result);
223 /*============================================================================*/
226 double ageAlongEdges(int node) { /* tree must be clock-like */
229 if(node>=0) return(0.);
231 for(i=0;i<tree_nedges;i++)
232 if(tree_lowerNodes[i]==node)
233 return(tree_edgeLengths[i]+ageAlongEdges(tree_upperNodes[i]));
238 /* compute set of parameters for a given clock-like tree */
239 void getExternalParams(double *result) {
242 for(i=0;i<tree_nedges;i++) {
243 if(tree_lowerNodes[i]<0&&tree_upperNodes[i]<0)
244 result[index[-tree_upperNodes[i]]]
245 =-log(ageAlongEdges(tree_lowerNodes[i])/ageAlongEdges(tree_upperNodes[i]));
250 /*============================================================================*/