src/nprsfunc.c

   1 /*
   2  *  nprsfunc.c
   3  *
   4  * (c) 2003  Gangolf Jobb and Korbinian Strimmer
   5  *
   6  *  Functions for nonparametric rate smoothing (NPRS)
   7  *  (see MJ Sanderson. 1997.  MBE 14:1218-1231)
   8  *
   9  *  This code may be distributed under the GNU GPL
  10  */
  11
  12
  13 #include <stdlib.h>
  14 #include <string.h>
  15 #include <math.h>
  16
  17 /*============================================================================*/
  18
  19 #define EPS 1e-6
  20 #define LARGENUM 1e99
  21
  22 /* original scale */
  23 #define MINPARAM EPS
  24 #define MAXPARAM 1-EPS
  25 #define STARTPARAM 0.5
  26
  27
  28 /* log scale */
  29 #define MINLPARAM log(MINPARAM)
  30 #define MAXLPARAM log(MAXPARAM)
  31 #define STARTLPARAM log(STARTPARAM)
  32
  33 #define ARRLEN 2048
  34 #define LABLEN 64
  35
  36 #define index aperates_index
  37
  38 int tree_lowerNodes[ARRLEN];
  39 int tree_upperNodes[ARRLEN];
  40 double tree_edgeLengths[ARRLEN];
  41 int tree_nedges;
  42 char tree_tipLabels[ARRLEN][LABLEN];
  43 int tree_ntips;
  44
  45 int nparams,index[ARRLEN],ancestor[ARRLEN];
  46
  47 /*============================================================================*/
  48
  49 void setTree(
  50  int *lowerNodes,
  51  int *upperNodes,
  52  double *edgeLengths,
  53  double *minEdgeLength,
  54  int *nedges,
  55  char **tipLabels,
  56  int *ntips,
  57  int *result
  58 ) {
  59  int i;
  60
  61  tree_nedges=*nedges;
  62
  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];
  67  }
  68
  69  tree_ntips=*ntips;
  70
  71  for(i=0;i<tree_ntips;i++) {
  72   strcpy(tree_tipLabels[i],tipLabels[i]);
  73  }
  74
  75  nparams=0;
  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];
  79  }
  80
  81  *result=0;
  82 }
  83
  84 /*============================================================================*/
  85
  86 void getNFreeParams(int *result) { *result=nparams; }
  87
  88 /*============================================================================*/
  89
  90 void getNEdges(int *result) { *result=tree_nedges; }
  91
  92 /*============================================================================*/
  93
  94 void getEdgeLengths(double *result) {
  95  int i;
  96
  97  for(i=0;i<tree_nedges;i++) result[i]=tree_edgeLengths[i];
  98
  99 }
 100
 101 /*============================================================================*/
 102
 103 double age(double *params,int node) {
 104  double prod;
 105
 106  if(node>=0) return(0.0);
 107  if(node==-1) return(1.0);
 108
 109  prod=0.0;
 110  while(node!=-1) {prod+=params[index[-node]]; node=ancestor[-node];}
 111
 112  return(exp(prod));
 113 }
 114
 115
 116 void getDurations(double *params,double *scale,double *result) {
 117  int i,low,upp;
 118
 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);
 123  }
 124
 125 }
 126
 127 /*============================================================================*/
 128
 129 void getRates(double *params,double *scale,double *result) {
 130  int i,low,upp;
 131
 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));
 136  }
 137
 138 }
 139
 140 /*============================================================================*/
 141
 142
 143 /* note on the choice of parameters:
 144  *
 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).
 149  *
 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
 152  * so we simply sum).
 153  */
 154
 155
 156 /* internal objective function - parameters are on log scale */
 157 void nprsObjectiveFunction(
 158  double *params,
 159  int *expo,
 160  double *result
 161 ) {
 162  int p,k,j;
 163  double me,rk,rj,sum,scale=1.0,durations[ARRLEN],rates[ARRLEN];
 164
 165  p=*expo;
 166
 167  getDurations(params,&scale,durations);
 168
 169  for(k=0;k<tree_nedges;k++) rates[k]=tree_edgeLengths[k]/durations[k];
 170
 171  me=0.; for(k=0;k<tree_nedges;k++) me+=rates[k]; me/=tree_nedges;
 172
 173  sum=0.;
 174
 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);
 179   }
 180  }
 181
 182  *result=sum;
 183
 184 }
 185
 186
 187 /* check parameter bounds on log scale */
 188 int checkLogParams(double *params)
 189 {
 190    int i;
 191    for(i=0; i<nparams; i++)
 192    {
 193      if(params[i] > MAXLPARAM || params[i] < MINLPARAM) return 0;  /* out of bounds */
 194    }
 195    return 1; /* within bounds */
 196 }
 197
 198
 199
 200 /*
 201  * public objective function
 202  * - parameters are on log scale
 203  * - if parameters are out of bounds function returns large value
 204  */
 205 void objFuncLogScale(
 206  double *params,
 207  int *expo,
 208  double *result
 209 )
 210 {
 211
 212  if( checkLogParams(params) == 0 ) /* out of bounds */
 213  {
 214    *result = LARGENUM;
 215  }
 216  else
 217  {
 218    nprsObjectiveFunction(params, expo, result);
 219  }
 220 }
 221
 222
 223 /*============================================================================*/
 224
 225
 226 double ageAlongEdges(int node) { /* tree must be clock-like */
 227  int i;
 228
 229  if(node>=0) return(0.);
 230
 231  for(i=0;i<tree_nedges;i++)
 232   if(tree_lowerNodes[i]==node)
 233    return(tree_edgeLengths[i]+ageAlongEdges(tree_upperNodes[i]));
 234
 235  return(0.);
 236 }
 237
 238 /* compute set of parameters for a given clock-like tree */
 239 void getExternalParams(double *result) {
 240  int i;
 241
 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]));
 246  }
 247
 248 }
 249
 250 /*============================================================================*/