src/delta_plot.c

   1 /* delta_plot.c       2011-06-23 */
   2
   3 /* Copyright 2010-2011 Emmanuel Paradis
   4
   5 /* This file is part of the R-package `ape'. */
   6 /* See the file ../COPYING for licensing issues. */
   7
   8 #include <R.h>
   9
  10 void delta_plot(double *D, int *size, int *nbins, int *counts, double *deltabar)
  11 {
  12         int x, y, u, v; /* same notation than in Holland et al. 2002 */
  13         int n = *size, nb = *nbins;
  14         double dxy, dxu, dxv, dyu, dyv, duv, A, B, C, delta;
  15         int xy, xu, xv, yu, yv, uv, i;
  16
  17         for (x = 0; x < n - 3; x++) {
  18                 xy = x*n - x*(x + 1)/2; /* do NOT factorize */
  19                 for (y = x + 1; y < n - 2; y++, xy++) {
  20                         yu = y*n - y*(y + 1)/2; /* do NOT factorize */
  21                         dxy = D[xy];
  22                         xu = xy + 1;
  23                         for (u = y + 1; u < n - 1; u++, xu++, yu++) {
  24                                 uv = u*n - u*(u + 1)/2; /* do NOT factorize */
  25                                 dxu = D[xu];
  26                                 dyu = D[yu];
  27                                 xv = xu + 1;
  28                                 yv = yu + 1;
  29                                 for (v = u + 1; v < n; v++, xv++, yv++, uv++) {
  30                                         dxv = D[xv];
  31                                         dyv = D[yv];
  32                                         duv = D[uv];
  33                                 /*      if (dxy <= dxu && dxy <= dxv && dxy <= dyu && dxy <= dyv && dxy <= duv) k=1; */
  34 /*                                      else if (duv <= dxy && duv <= dxu && duv <= dxv && duv <= dyu && duv <= dyv) k=1; */
  35 /*                                      else if (dxu <= dxy && dxu <= dxv && dxu <= dyu && dxu <= dyv && dxu <= duv) k=2; */
  36 /*                                      else if (dyv <= dxy && dyv <= dxu && dyv <= dxv && dyv <= dyu && dyv <= duv) k=2; */
  37 /*                                      else if (dxv <= dxy && dxv <= dxu && dxv <= dyu && dxv <= dyv && dxv <= duv) k=3; */
  38 /*                                      else if (dyu <= dxy && dyu <= dxu && dyu <= dxv && dyu <= dyv && dyu <= duv) k=3; */
  39                                         //Rprintf("%d\t%d\t%d\t%d\t%d\t%d\t\n", xy, xu, xv, yu, yv, uv);
  40                                         //Rprintf("dxy = %f\tdxu = %f\tdyu = %f\n", dxy, dxu, dyu);
  41                                         //Rprintf("D[xv] = %f\tD[yv] = %f\tD[uv] = %f\n", D[xv], D[yv], D[uv]);
  42                                         /* switch (k) { */
  43 /*                                      case 1 : A = dxv + dyu; B = dxu + dyv; C = dxy + duv; break; */
  44 /*                                      case 2 : A = dxv + dyu; B = dxy + duv; C = dxu + dyv; break; */
  45 /*                                      case 3 : A = dxu + dyv; B = dxy + duv; C = dxv + dyu; break; */
  46 /*                                      } */
  47                                         //Rprintf("A = %f\tB = %f\tC = %f\n", A, B, C);
  48                                         A = dxv + dyu; B = dxu + dyv; C = dxy + duv;
  49                                         //Rprintf("A = %f\tB = %f\tC = %f\n", A, B, C);
  50                                         if (A == B && B == C) delta = 0; else while (1) {
  51                                                 if (C <= B && B <= A) {delta = (A - B)/(A - C); break;}
  52                                                 if (B <= C && C <= A) {delta = (A - C)/(A - B); break;}
  53                                                 if (A <= C && C <= B) {delta = (B - C)/(B - A); break;}
  54                                                 if (C <= A && A <= B) {delta = (B - A)/(B - C); break;}
  55                                                 if (A <= B && B <= C) {delta = (C - B)/(C - A); break;}
  56                                                 if (B <= A && A <= C) {delta = (C - A)/(C - B); break;}
  57                                         }
  58                                         /* if (delta == 1) i = nb - 1; else */
  59                                         i = delta * nb;
  60                                         //Rprintf("delta = %f\ti = %d\n", delta, i);
  61                                         counts[i] += 1;
  62                                         deltabar[x] += delta;
  63                                         deltabar[y] += delta;
  64                                         deltabar[u] += delta;
  65                                         deltabar[v] += delta;
  66                                 }
  67                         }
  68                 }
  69         }
  70 }