# Purpose: Produce Nucleotide-Nucleotide match from NPMatches
# Author:  Lukas Knecht
# Created: 25 Apr 1995
#
NNMatch := proc (m1: NucPepMatch, m2: NucPepMatch, df1: database, df2: database)
  global NucDB;
  description 'Calculates the Nucleotide-Nucleotide alignment implied by
  m1 and m2. m1 and m2 must refer to the same protein.';

  ForwardPep := proc (spm: [string, string], f: integer, t: integer)
    o := f; j := 0;
    for i while o < t do
      if spm[1,i] = '_' then o := o + 1
      elif spm[2,i] <> '_' then
	j := j + 1; if j = 3 then j := 0; o := o + 1 fi
      fi
    od;
    i
  end;

  oldNDF := NucDB;
  if nargs > 2 then NucDB := df1 fi;
  spm1 := DynProgNucPepString (m1);
  if nargs > 3 then NucDB := df2 fi;
  spm2 := DynProgNucPepString (m2);
  i1 := ForwardPep (spm1, m1[PepOffset], m2[PepOffset]);
  i2 := ForwardPep (spm2, m2[PepOffset], m1[PepOffset]);
  res1 := CreateString (length (spm1[1]) + length (spm2[1])); res2 := copy (res1);
  p1 := CreateArray (1..3); p2 := copy (p1);
  l := 0;
  for o from max (m1[PepOffset], m2[PepOffset]) to
  min (m1[PepOffset] + m1[PepLength], m2[PepOffset] + m2[PepLength]) - 1 do
    # Check for amino acid deletions
    if spm1[1,i1] = '_' then
      i1 := i1 + 1;
      if spm2[1,i2] = '_' then
	i2 := i2 + 1
      else
	j := 0;
	for i2 from i2 while j < 3 do
	  l := l + 1; res1[l] := '_'; res2[l] := spm2[1,i2];
	  if spm2[2,i2] <> '_' then j := j + 1 fi
	od
      fi
    elif spm2[1,i2] = '_' then
      i2 := i2 + 1;
      j := 0;
      for i1 from i1 while j < 3 do
	l := l + 1; res1[l] := spm1[1,i1]; res2[l] := '_';
	if spm1[2,i2] <> '_' then j := j + 1 fi
      od
    else
      # Get matching bases from m1
      j := 0;
      for i from i1 while j < 3 do
	if spm1[2,i] <> '_' then j := j + 1; p1[j] := i fi
      od;
      # Get matching bases from m2
      j := 0;
      for i from i2 while j < 3 do
	if spm2[2,i] <> '_' then j := j + 1; p2[j] := i fi
      od;
      # Encode strings
      for j to 3 do
	for i1 from i1 to p1[j] - 1 do
	  l := l + 1; res1[l] := spm1[1,i1]; res2[l] := '_'
	od;
	for i2 from i2 to p2[j] - 1 do
	  l := l + 1; res1[l] := '_'; res2[l] := spm2[1,i2]
	od;
	l := l + 1; res1[l] := spm1[1,i1]; res2[l] := spm2[1,i2];
	i1 := i1 + 1; i2 := i2 + 1
      od
    fi
  od;
  NucDB := oldNDF;
  [res1[1..l], res2[1..l]]
end: