cpnni.pas

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

      LeafCount := LeafCount + S1^.Treeleaves;

      {2:  Red subtree }
      new (S2, Init);
      S2^.Root := There^.Child^.Sib;
      S2^.Root^.anc := NIL;
      S2^.Root^.sib := NIL;
      S2^.SetLeaves (S2^.Root^.GetWeight);
      S2^.SetInternals(Pred (S2^.TreeLeaves));

      LeafCount := LeafCount + S2^.Treeleaves;

      { Ensure green subtree's root doesn't still
        have red subtree as a sibling }
      S1^.Root^.sib := NIL;

      AddExtraNodeGR (S1^);  { do this here }
      AddExtraNodeGR (S2^);


      { 3: Blue subtree }
      { Make There the extra node }
      with There^ do begin
         SetLeaf;
         SetNamePtr (ExtraNodeCount);
         Child := NIL;
         end;
      T.SetLeaves (T.TreeLeaves - LeafCount + 1);
      T.SetInternals (Pred (T.TreeLeaves));
      T.NL[ExtraNodeCount] := There;
      end;

   { Fix up sub trees

     The resulting subtree resembles a pruned tree in that
     the NL fields are indexed by the order of the values
     in the original tree, e.g.


     1  2  3  4  5  6  NL
     ---------------
     -  -  3  -  5  6  NamePtr
     -  -  3  -  1  2  LeafNum

     5   6   3
      \   \ /
       \   * [2,3]
        \ /
         * [1..3]

     Hence, clusters in tree reflect number of leaves in subtree,
     but NamePtr fields reflect original profile. This is essential
     to ensure that leaves are consistently labeled.
   }


   T.ReRoot (T.NL[ExtraNodeCount]);
   new (S0, Init);
   S0^.Root := T.CopyOfSubTree (T.Root);
   S0^.SetLeaves (T.TreeLeaves);
   S0^.SetINternals (T.TreeInternals);
   S0^.BuildLeafList;
   {$IFDEF NNIDEbug}
   writeln (output, 'Ordering S0');
   {$ENDIF}

   S0^.OrderLeaves;
   S0^.BuildClusters;
   S0^.CorrectWeights;

   S1^.BuildLeafList;
   {$IFDEF NNIDEbug}
   writeln (output, 'Ordering S1');
   {$ENDIF}

   S1^.OrderLeaves;
   S1^.BuildClusters;
   S1^.CorrectWeights;
   if not Edge then begin
      S2^.BuildLeafList;
      {$IFDEF NNIDEbug}
      writeln (output, 'Ordering S2');
      {$ENDIF}

      S2^.OrderLeaves;
      S2^.BuildClusters;
      S2^.CorrectWeights;
      end;

   { Store in data structure }
   S[0] := S0;
   S[1] := S1;
   if not Edge then
      S[2] := S2;
end;

{-----------------------------MarkASubTree---------------------------------}

{ For each leaf in T not of the desired Color,
  mark it and its first nonmarked ancestor.
  This is logically equivalent to pruning the
  tree of all leaves except those "of Color." }
procedure MarkASubTree (var T:TREEOBJ;Color:Byte);
var
   i, j: integer;
   q : NODE_PTR;
begin
   { Remember that the leaves of a tree need not be
     labeled 1..n because of pruning. }
   j := 0;
   i := 0;
   while (j < T.TreeLeaves) do begin
      repeat
         Inc (i);
      until (T.NL[i] <> NIL);
      if (T.NL[i]^.NodeIndex <> Color) then begin
         T.NL[i]^.Mark;
         q := T.NL[i]^.NodeAncestor;
         while q^.IsMarked do
            q := q^.NodeAncestor;
         q^.Mark;
         end;
      Inc (j);
      end;
end;

{-----------------------------GetSubTree-----------------------------------}

{ Given a start tree T return in SubT a subtree. Nodes belonging in
  the subtree are unmarked, all other nodes are marked. The subtree
  may be imbedded in T, e.g.:

  T =

  B   B   G   B   G   R G   R

      +---LimitSister
      |
  5   6   4   3   1   2 7   8
   \   \   .   \   . /   . /
    \   \   .   \   *     *
     \   \   .   \ .     .
      \   \   .   *     .
       \   \   .   .   .
        \   \   .   . .
         \   \   .   o
          \   \   . .
           \   \   o <-- (Limit)
            \   \ /
             \   *
              \ /
               *

  o = unmarked, . = imbedded subtree.

  The Green subtree returned is

  SubT =

  x   4   1   7
   .   \   \ /
    .   \   *
     .   \ /
      .   *
       . /
        *

  Algorithm operates by traversing T and using SubT's stack to
  keep track of nodes added to the subtree. The first node
  added may be just a place holder below SubT's root,
  and hence may be pruned off the subtree.


  The original T IS NOT altered or destroyed by this procedure,
  but all nodes are unmarked.

}
procedure GetSubTree (var T, SubT: TREEOBJ; Limit: NODE_PTR; Color: Byte);
var
   BlueOGPtr : NODE_PTR; { first node in blue subtree above Limit }
   GetBlueOG : Boolean;  { flag to trap BlueOGPtr }
   BlueSister: NODE_PTR;
   LimitSister : NODE_PTR;

   procedure Traverse (p:NODE_PTR);
   var
      q, r: NODE_PTR;
   begin
      if (p <> NIL) then begin
         if (p = Limit) then
            { We've passed Limit so prepare to trap the first
              node that is in a subtree...}
            GetBlueOG := True;
         if not p^.IsMarked then begin
            { p is part of subtree }
            q := SubT.ST.TopOfStack;
            if (q^.Child = NIL) then begin
               r := new (NODE_PTR, Init);
               q^.child := r;
               r^.anc := q;
               if (p^.IsLeaf) then begin
                  SubT.IncLeaves;
                  r^.SetLeaf;
                  r^.SetNamePtr (p^.NodeNamePtr);
                  SubT.NL[p^.NodeNamePtr] := r;
                  end
               else begin
                  r^.SetWeight (p^.GetWeight);
                  SubT.IncInternals;
                  SubT.ST.Push (r);
                  end;
               end
            else begin
               r := new (NODE_PTR, Init);
               q^.child^.sib := r;
               r^.anc := q;
               if (p^.IsLeaf) then begin
                  SubT.IncLeaves;
                  r^.SetLeaf;
                  r^.SetNamePtr (p^.NodeNamePtr);
                  SubT.NL[p^.NodeNamePtr] := r;
                  end
               else begin
                  r^.SetWeight (p^.GetWeight);
                  SubT.IncInternals;
                  SubT.ST.Push (r);
                  end;
               end;
            { ...if trap is set, capture the first node
              in the subtree that occurs above Limit and
              release the trap. }
            if GetBlueOG then begin
               BlueOGPtr := r;
               GetBlueOG := false;
               end;
            { trap the first blue node in SubT that corresponds to
              the sister of limit in T }
            if (Color = B) and (p = LimitSister) then
               BlueSister := r;
            end;

         Traverse (p^.child);
         if (not p^.IsLeaf) and (not p^.IsMarked) then
            SubT.ST.Pop;
         { Unmark as we go down tree }
         p^.UnMark;
         { Release trap so that any nodes
           BELOW limit aren't spuriously
           caught. }
         if (p = Limit) then
            GetBlueOG := false;
         Traverse (p^.sib);
         end;
   end;

begin
   BlueOGPtr  := NIL;
   GetBlueOG  := false;

   { Get the sister of limit on T }
   if Limit = Limit^.anc^.child then
      LimitSister := Limit^.sib
   else LimitSister := Limit^.anc^.child;

   SubT.Init;
   SubT.MakeRoot;
   SubT.IncInternals;
   SubT.ST.Push (SubT.Root);
   Traverse (T.Root);

   if (Color <> B) then begin
      { Add extra node }
      SubT.Root^.child^.sib := ExtraNode;
      SubT.Root^.child^.sib^.anc := SubT.Root;
      SubT.IncLeaves;
      SubT.NL[SubT.Root^.child^.sib^.NodeNamePtr] := SubT.Root^.child^.sib;
      end
   else begin
      { prune off subroot }
      SubT.Root := SubT.Root^.child;
      Dispose (SubT.Root^.anc, Done);
      SubT.Root^.anc := NIL;

      { Insert an extra node }
      if (BlueOGPtr <> NIL) then begin
         { >= 1 blue leaves occur above Limit, add the
           Extra node below the first one. }
         AddExtraNodeB (SubT, BlueOGPtr);
         end
      else begin
         { Blue leaves are all below Limit, so
           add extra node below Limit's sister
           node (which must be blue). }
         AddExtraNodeB (SubT, BlueSister);
         end;

      { Ensure consistent rooting by using extra }
      SubT.ReRoot (SubT.NL[ExtraNodeCount]);

      end;

   { Fix up sub tree

     The resulting subtree resembles a pruned tree in that
     the NL fields are indexed by the order of the values
     in the original tree, e.g.


     1  2  3  4  5  6  NL
     ---------------
     -  -  3  -  5  6  NamePtr
     -  -  3  -  1  2  LeafNum

     5   6   3
      \   \ /
       \   * [2,3]
        \ /
         * [1..3]

     Hence, clusters in tree reflect number of leaves in subtree,
     but NamePtr fields refelct original profile. This is essential
     to ensure that leaves are consistently labeled.
   }


{$IFDEF NNIDEbug}
   writeln (output, 'SubT.OrderLeaves');
{$ENDIF}

   SubT.OrderLeaves;      { Ensure LeafNum field allocated }
   SubT.BuildClusters;    { Ensure valid clusters }
   SubT.CorrectWeights;   { Ensure correct weights }
end;

{-----------------------------DRA------------------------------------------}

{ Recursively compute NNI between T1 and T2 }

procedure doDRA (var T1, T2: TREEOBJ);
var
   S1, S2: SUBTREES;
   i, c, k : integer;
   There, Limit : NODE_PTR;
   Edge : Boolean;
   m1, m2: Longint;
begin
   if (T1.TreeLeaves < 4) then
      Exit;

   { Color the goal tree }
   ColorTree (T2, There, Edge);
   if Edge then
      k := 1
   else k := 2;

   { Color the start tree }
   LabelT (T1, T2);

   { Return the NNI needed to transform start- into goal-tree }
   c := Count (T1, Limit);
   NNICount := NNICount + c;

   { Decompose trees... }
   { Goal tree }

   {$IFDEF NNIDEbug}
   writeln (output, 'Doing goal subtrees');
   {$ENDIF}

   DecompGoalTree (T2, S2, There, Edge);

   {$IFDEF NNIDEbug}
   writeln (output, 'Goal subtrees');
   S2[0]^.WriteTree (output);
   S2[1]^.WriteTree (output);
   if not Edge then
      S2[2]^.WriteTree (output);
   {$ENDIF}

   { Start tree }
   { blue }
   {$IFDEF NNIDebug}
   writeln (output, 'Blue');
   {$ENDIF}

   MarkASubTree (T1, 0);
   new (S1[0]);
   GetSubTree (T1, S1[0]^, Limit, 0);

   { green }
   {$IFDEF NNIDebug}
   writeln (output, 'Green');
   {$ENDIF}

   MarkASubTree (T1, 1);
   new (S1[1]);
   GetSubTree (T1, S1[1]^, Limit, 1);

   { red }
   if not Edge then begin
      {$IFDEF NNIDebug}
      writeln (output, 'Red');
      {$ENDIF}
      MarkASubTree (T1, 2);
      new (S1[2]);
      GetSubTree (T1, S1[2]^, Limit, 2);
      end;

   {$IFDEF NNIDebug}
   writeln (output, 'Start subtrees');
   S1[0]^.WriteTree (output);
   S1[1]^.WriteTree (output);
   if not Edge then
      S1[2]^.WriteTree (output);
   {$ENDIF}


   { recursive loop }
   Inc (ExtraNodeCount);
   Inc (LoopCount);
   for i := 0 to k do begin
      doDRA (S1[i]^, S2[i]^);

      { Free the start subtrees }
      Dispose (S1[i], Done);
      Dispose (S2[i], Done);
      end;
   Dec (ExtraNodeCount);
   Dec (LoopCount);
end;

{-----------------------------NNI------------------------------------------}

{ Return dra between T1 and T2, destroying T1 and T2
  in the process. }
function DRA (var T1, T2: TREEOBJ; rt:integer):integer;
var
   m1, m2: longint;
begin
   {$IFDEF NNIDEbug}
   writeln (output, '*****In NNI*****');
   {$ENDIF}


   NNICount  := 0;
   LoopCount := 0;

   { ExtraNodeCount must start at the first cell in
     both T1 and T2 NL array beyond the leaves in
     those trees. Need to take into account that
     (1) T1 and T2 may be pruned, and
     (2) if they're from different profiles then
         they may be unequally pruned.
   }
   ExtraNodeCount := Max (Succ(T1.TreeLeaves + T1.TreePrunedleaves),
                          Succ(T2.TreeLeaves + T2.TreePrunedleaves));

   { Root T1 and T2 with same leaf }
   T1.ReRoot (T1.Ithleaf (rt));
   T2.ReRoot (T2.IthLeaf (rt));

   { Compute }
   doDRA (t1, T2);

   { Dispose of the remains of T1 and T2 }
   T1.Done;
   T2.Done;

   { Return }
   DRA := NNICount;

   {$IFDEF NNIDEbug}
   writeln (output, '*****Out of NNI*****');
   {$ENDIF}

end;

begin
end.