cprtree.pas

生物信息学中的遗传数据分析的delphi源码。

{*********************************************}
{                                             }
{    COMPONENT for MS DOS and MS WINDOWS      }
{                                             }
{    Source code for Turbo Pascal 6.0 and     }
{    Turbo Pasacal for Windows 1.0 compilers. }
{                                             }
{    (c) 1991, Roderic D. M. Page             }
{                                             }
{*********************************************}


{$I CPDIR.INC}

unit cprtree;

{
   Random tree generators

   1. Uniform unlabeled binary rooted trees
   2.    "        "       "    unrooted trees
   3. All rooted/unrooted, binary, labeled trees

  The methods for unlabeled trees makes extensive use of Furnas' (1984)
  paper. The routines here and in CPTREE.PAS are his, and also allow
  for unique labelling of both rooted and unrooted tree topologies
  using the RRANK and CRANK functions in CPTREE.

  26 Nov 1991 Sketched out.

*}

interface

uses
   cpheader, { Resource ids }
   cpwcdial,
   cpwvars,
   cperror,
   cputil,   { utilities }
   cptree,   { tree }
   cpanc,
   cpbtree,
   cplabels,
   cpwntbuf;

const
   rt_Rooted      = $0001;
   rt_Labeled     = $0002;
   rt_Dendrogram  = $0004;
   rt_Exhaustive  = $0008;
   rt_Sample      = $0010;
   rt_Binary      = $0020;

   rt_HDU = rt_Labeled or rt_Binary;
   rt_HDR = rt_HDU or rt_Rooted;      { = HD }
   rt_HEU = rt_Binary;
   rt_HER = rt_HEU or rt_Rooted;      { = HE }
   rt_HM  = rt_HDR or rt_Dendrogram;  { = HM }

type
   ALLNTREES_PTR = ^ALLNTREES;
   { Pointer to [ALLNTREES] }
   ALLNTREES = object (ALLTREES)
   { Object to generate all rooted and unrooted, labelel n-trees }
      B : PTREE_BUFFEROBJ; {Pointer to tree buffer to store trees }
      constructor Init (ABuffer: PTREE_BUFFEROBJ;L:LABELOBJ_PTR;
                        leaves:integer);
         {Call [ALLTREES], set [B]
          to <\b ABuffer>, and Count (private) to 0.}
      procedure TreeOp;virtual;
         {Store current tree in buffer and update counter display}
      private
      Count:integer;
      end;



procedure UniformUnlabeledRooted (rank, leaves:integer; var T:TREEOBJ);
{Return in <\b T> the unlabeled rooted tree of size <\b leaves> with the
  <\b rank>th topology in the list of ordered topologies. } 

procedure UniformUnlabeledUnrooted (rank, leaves:integer; var T:TREEOBJ;
                                    var phin: longint);
{Return in <\b T> the unlabeled unrooted tree of size <\b leaves> with the
 <\b rank>th topology in the list of ordered unrooted topologies. Must pass
 the value of [phi](n) in <\b phin>} 

implementation

type
   func = function (i,z,Vab:real):real;

var
   ternod,
   intnod,
   root : integer;
   AF   : ANCFUNCOBJ;



{-----------------------------f1-------------------------------------------}

{ Furnas (1984:221) equation for i1 }
function f1 (i,z,Vab:real):real;far;
begin
   f1 := Sqr(i)*i - 3 * (z+1)*Sqr(i) + (3*sqr(z)+6*z+2)*i - 6*Vab;
end;

{-----------------------------f2-------------------------------------------}

{ Furnas (1984:221) equation for i2 }
function f2 (i,z,Vab:real):real;far;
begin
   f2 := Sqr(i)*i - 3*z*Sqr(i) + (3*sqr(z)-1)*i + 3*z*(z+1)- 6*(Vab+1);
end;

{-----------------------------RootOf---------------------------------------}

{ Find the root of function Operation. Used to solve the cubic
  formulas f1 and f2 above (see Furnas, 1984:221).

  Code taken from T. Rugg & P. Feldman. 1986.
  "Turbo Pascal program library," Que, Indianapolis.
}
procedure RootOf (Operation: Func; LowX, HighX, z, Vab:real;
                  var Result:real; var OK:Boolean);
const
   MaxNumSegs = 500;
   MaxError   = 1.0E-06;
var
   X, Xlo, Xhi, Flo, Fhi, Width, Test : real;
   j, numsegs : integer;

   procedure Iterate;
   begin
      repeat
         X := (FHi * XLo - Flo * XHi) / (Fhi-Flo);
         Test := Operation (X, z,vab);
         if (abs (Test) < MaxError) then begin
            Result := X;
            OK     := true;
            end
         else
            if (Test * flo) > 0.0 then begin
               xlo := x;
               flo := test;
               end
            else begin
               xhi := x;
               fhi := test;
               end
      until OK;
   end;

begin
   OK := False;
   xlo := LowX;
   xhi := HighX;
   flo := Operation (xlo, z, vab);
   fhi := Operation (xhi, z, vab);
   if (abs(flo) < MaxError) then begin
      OK := true;
      Result := xlo;
      exit;
      end;
   if (abs(fhi) < MaxError) then begin
      OK := true;
      Result := xhi;
      exit;
      end;
   if (flo*fhi) < 0.0 then begin
      iterate;
      exit;
      end;
   NumSegs := 2;
   repeat
      width := (xhi-xlo)/NumSegs;
      for j := 1 to (NumSegs div 2) do begin
         x := xlo + width * (2 * j - 1);
         fhi := Operation (x, z, vab);
         if (fhi*flo) < 0.0 then begin
            xhi := x;
            iterate;
            exit;
            end;
         end;
      numsegs := numsegs * 2;
   until (Numsegs > maxNumSegs);
end;

{-----------------------------iW1------------------------------------------}

{ Furnas' (1984:214, eqn 9) formula for solving W(x,i,j) for i }
function iW1 (x,w:longint):longint;
begin
   iW1 := max (Trunc ((2*x-1-Sqrt(Sqr(2*x-1) - 8*(w-x+1)))/2),
               Trunc ((2*x+1-Sqrt(Sqr(2*x+1) - 8*w))/2));
end;

{-----------------------------iW2------------------------------------------}

{ Furnas' (1984:214, eqn 10) formula for solving W(x,i,j) for j }
function iW2 (x,w,w1:longint):longint;
begin
   iW2 := w - d(x) + d(x-W1) + W1;
end;

{-----------------------------invrank--------------------------------------}

{ Furnas' (1984:215,231) recursive procedure to generate
  a tree with the LLR toplogy given by rank.

  The first call requires that
     ternod=0;
     root=intnod=n+1;
     AF.Init;

  The result is an ancestor function AF containing a
  tree with the desired topology. The leaves are always
  in 1,..,n order, so a random tree will require relabeling
  of the leaves.
}

procedure invrank (rank:longint; n:integer; node:integer);
var
   a,               { size of left and right subtrees }
   b     : integer;
   bl,              { no. of trees preceding left and right subtrees }
   br: longint;
   x,
   p,
   V     : longint; { local variables }
   lnode,           { indices of left and right subtrees }
   rnode : integer;

begin
   if (n > 1) then begin
      { Find size of left and right subtrees }
      a := 0;
      p := 0;
      repeat
         Inc (a);
         p := p + Z[a]*Z[n-a];
      until (p >= rank);
      p := p - Z[a]*Z[n-a];
      b := n - a;

      { compute relative ranking of subtrees }
      V := rank - p - 1;
      if (a = b) then begin
         x  := Z[n div 2];
         bl := iW1 (x,V);
         br := iW2 (x,V,bl);
         end
      else begin
         bl := V div Z[b];
         br := V mod Z[b];
         end;

      { create subtrees }
      { make left desc }
      if (a > 1) then begin
         { internal }
         Inc (intnod);
         lnode := intnod;
         end
      else begin
         { terminal }
         Inc (ternod);
         lnode := ternod;
         end;

      { place in ancestor function }
      AF.A[lnode+2] := node;

      { generate left subtree }
      invrank (bl+1, a, lnode);

      { make right desc }
      if (b > 1) then begin
         { internal }
         Inc (intnod);
         rnode := intnod;
         end
      else begin
         { terminal }
         Inc (ternod);
         rnode := ternod;
         end;

      { place in ancestor function }
      AF.A[rnode+2] := node;

      { generate right subtree }
      invrank (br+1, b, rnode);
      end;
end; {invrank}

{-----------------------------UniformUnlabeledRooted-----------------------}

{ Return in T an "unlabeled" tree with the rankth topology for
  rooted binary trees with a given number of leaves.

  Although unlabeled, the leaves are ordered left to right
  <1,..,leaves>. Rank is in the range 1..Z[leaves].
}
procedure UniformUnlabeledRooted (rank, leaves:integer; var T:TREEOBJ);
begin
   ternod := 0;
   intnod := leaves + 1;
   root   := intnod;
   AF.Init;
   AF.A[1] := leaves;
   AF.A[2] := Pred(Leaves);
   invrank (rank, leaves, root);
   AF.AncFuncToTree (T);
end;


{-----------------------------UniformUnlabeledUnrooted---------------------}

{ Return in T an "unlabeled" tree with the rankth topology for
  unrooted binary trees with a given number of leaves.

  Although unlabeled, the leaves are ordered left to right
  <1,..,leaves>. Rank is in the range 1..Z[leaves].

  Routine requires that the number of unrooted triple trunk
  tree topologies be already calculated (phin). This is
  done by calling phi(leaves) before any calls to this
  proc. (This is simply to avoid repetetive calls to phi(leaves).
}

procedure UniformUnlabeledUnrooted (rank, leaves:integer; var T:TREEOBJ;
                                    var phin: longint);
var
   a,           { sizes of subtrees }
   b,
   c : integer;
   bl,
   bm,
   br : longint; { no. of trees before @subtree in LLR order }
   Vab,
   w,
   x1 : longint;
   i1,
   i2 : real;
   OK :Boolean;
begin
   if (rank <= phin) then begin
      { triple trunk tree }
      { i. get subtree sizes }
      a  := 0;
      x1 := 0;
      while (x1 < rank) do begin
         Inc (a);
         x1 := U(leaves,a,a);
         end;
      a := a - 1;

      b  := 0;
      x1 := 0;
      while (x1 < rank) do begin
         Inc (b);
         x1 := U(leaves,a,b);
         end;
      b := b - 1;
      c := leaves - a - b;

      { ii. ranks of subtrees }
      Vab := rank - U(leaves,a,b) - 1;
      if (a<b) then begin
         if (b<c) then begin                         { a<b<c}
            bl := Trunc (Vab/(Z[b]*Z[c]));
            bm := Trunc ((Vab-bl*Z[b]*Z[c])/Z[c]);
            br := Vab mod Z[c];
            end
         else begin                                  {a<b=c}
            bl := Trunc (Vab/d(Z[b]));
            w  := Vab-bl*d(Z[b]);
            bm := iW1 (Z[b],w);
            br := iW2 (Z[b], w, bm);
            end;
         end
      else begin
         if (b<c) then begin                         { a=b<c}
            w  := Trunc (Vab/Z[c]);
            bl := iW1 (Z[a], w);
            bm := iW2 (Z[a], w, bl);
            br := Vab mod Z[c];
            end
         else begin                                  {a=b=c}
            Rootof (f1, 0 ,leaves div 2, Z[a], Vab, i1, OK);
            RootOf (f2, 0 ,leaves div 2, Z[a], Vab, i2, OK);
            bl := max (Trunc(i1),Trunc(i2));
            w  := Vab-d3(Z[a])+d3(Z[a]-bl);
            bm := bl + iW1(Z[a]-bl,w);
            br := bl + iW2(Z[a]-bl,w,bm-bl);
            end;
         end;
      end
   else begin
      { double trunk tree }
      { i. get subtree sizes }
      a := leaves div 2;
      c := leaves - a;

      { ii. ranks of subtrees }
      w := rank - phin - 1;
      bl := iW1 (Z[a], w);
      br := iW2 (Z[a], w, bl);
      bm := -1;
      end;

   { build tree }

   ternod  := 0;
   intnod  := leaves + 1;
   root    := intnod;
   AF.Init;
   AF.A[1] := leaves;

   { left subtree }
   if (a = 1) then begin
      { leaf }
      Inc (ternod);
      AF.A[ternod + 2] := root;
      end
   else begin
      { internal }
      Inc (intnod);
      AF.A[intnod + 2] := root;
      invrank (bl + 1, a, intnod);
      end;
   if (bm > -1) then begin
      { middle subtree }
      if (b=1) then begin
         { leaf }
         Inc (ternod);
         AF.A[ternod + 2] := root;
         end
      else begin
         { internal }
         Inc(intnod);
         AF.A[intnod + 2] := root;
         invrank (bm + 1, b, intnod);
         end;
      AF.A[2] := leaves-2; { triple trunk tree }
      end
   else
      AF.A[2] := leaves-1; { double trunk tree }

   { right subtree }
   if (c=1) then begin
      { leaf }
      Inc (ternod);
      AF.A[ternod + 2] := root;
      end
   else begin
      { internal }
      Inc (intnod);
      AF.A[intnod + 2] := root;
      invrank (br + 1, c, intnod);
      end;

  AF.AncFuncToTree (T);
  if (bm > -1) then
     T.MakeRootBinary;
end;



   constructor ALLNTREES.Init(ABuffer : PTREE_BUFFEROBJ;
                              L:LABELOBJ_PTR;leaves:integer);
   begin
      B := ABuffer;
      ALLTREES.Init (leaves, L);
      Count := 0;
   end;

   procedure ALLNTREES.TreeOp;
   begin
      Inc (Count);
      { Store T in P's buffer }
      B^.PutTree (T);
      Error := B^.BufferError;

      {$IFDEF WINDOWS}
      if (Error=erOK) and (Counter <> NIL) then begin
         Counter^.UpDateMeter (id_Meter, Count);
         Counter^.UpDateNumber (id_TreesRead, Count);
         if bUserAbort then
            Error := erUserAbort;
         end;
      {$ENDIF}
   end;




end.