#! /usr/bin/gawk -f
# Last edited on 2005-10-17 17:13:59 by stolfi

BEGIN {
  usage = ( \
    "simula-dna \\\n" \
    "   -v tab1=ARQTAB1.frq \\\n" \
    "   -v tab3=ARQTAB3.frq \\\n" \
    "   -v nbases=NUM  \\\n" \
    "   -v tmin=NUM -v tc=NUM -v tn=NUM \\\n" \
    "   -v out=NAME" \
  );
  # Simula a geração de DNA segundo nosso modelo.
  # O gerador de seqencias não codificadoras (rótulo 'N') simplesmente 
  # concatena bases independentemente, com distribuição da tabela {tab1}.
  # O gerador de seqüências codificadoras (rótulos 'D', 'E', 'F')]
  # contatena triplas de bases com a distribuição da tabela {tab3}.
  # A saída dos dois geradores é recortada em pedaços de tamanho 
  # míimo {tmin} e tamanho médio {tn} e {tc}, respectivamente, e 
  # esses pedaços são concatenados alternadamente.

  if (tab1 == "") { arg_error("deve definir tab1"); }
  if (tab3 == "") { arg_error("deve definir tab3"); }
  if (nbases == "") { arg_error("deve definir nbases"); }
  if (tmin == "") { arg_error("deve definir tmin"); }
  if (tc == "") { arg_error("deve definir tc"); }
  if (tn == "") { arg_error("deve definir tn"); }
  if (out == "") { arg_error("deve definir out"); }

  arqbas = ( out ".bas" );  # Arquivo de bases
  arqrot = ( out ".rot" );  # Arquivo de rotarito

  split("", frq1); # Define {frq1} como matriz.
  le_tabela(tab1, frq1, 0);

  split("", frq3);
  le_tabela(tab3, frq3, 1);

  # Gera uma seqüência não codificadora com pelo menos {2*nbases} bases:
  non = "";
  for (i = 0; i < 2*nbases; i++)
    { c = gera_cadeia(frq1); non = ( non c ); }

  # gera uma seqüência codificadora com pelo menos {2*nbases} bases:
  cod = "";
  for (i = 0; i < 2*nbases; i += 3)
    { c = gera_cadeia(frq3); cod = ( cod c ); }

  # recorta e junta:
  seq = "";
  rot = "";
  fase = 0;
  while (length(seq) < nbases)
    { 
      # Pega um pedaço da seqüência {non}:
      m = gera_comprimento(tmin,tn,nbases);
      ped = substr(non, 1, m); non = substr(non, m+1);
      seq = ( seq ped );
      for(i = 0; i < m; i++) { rot = (rot "I"); }

      # Pega um pedaço da seqüência {cod}:
      m = gera_comprimento(tmin,tc,nbases);
      ped = substr(cod, 1, m); cod = substr(cod, m+1);
      seq = ( seq ped );
      for(i = 0; i < m; i++) { rot = (rot substr("DEF", fase+1, 1)); fase = (fase + 1) % 3; }
    }

  # Escreve as primeiras {nbases} letras de {seq} e {rot}:
  for (i = 0; i < nbases; i++)
    { 
      if ((i > 0) && (i % 60 == 0)) 
        { printf "\n" > arqbas;  printf "\n" > arqrot; }
      printf "%s", substr(seq, i+1, 1) > arqbas;
      printf "%s", substr(rot, i+1, 1) > arqrot;
    }
  printf "\n" > arqbas;  close(arqbas);
  printf "\n" > arqrot;  close(arqrot);
        
}

function le_tabela(nome,frq,codif,     lin,fld,nfld,ok,tot,s,nok)
{
  # Le uma tabela de freqüências do arquivo {nome}, coloca em {frq[s]} a
  # probabilidade da tupla {s}.  Se {codif=1}, considera apenas tuplas
  # com evento "DEF"; se {codif=0}, considera apenas aquelas com evento "I".
  printf "lendo tabela \"%s\"...\n", nome > "/dev/stderr";
  tot = 0;
  ERRNO = "";
  while ( (getline lin < nome) > 0)
    {
      if (lin !~ /^[>]/)
        {
           nfld = split(lin, fld);
           if (nfld != 4) { table_error(nome, ("bad num fields")); }
           if (codif)
             { ok = (fld[1] == "EFG"); }
           else 
             { ok = (fld[1] == "I"); }
           if (ok)
             { frq[fld[2]] = fld[3]; tot += fld[3]; nok++; }
        }
    }
  if (ERRNO != "") { table_error(nome, ("erro de leitura ERRNO = \"" ERRNO "\"")); }
  if (tot == 0) { table_error(nome, ("frequencia total nula")); }
  # Converte contagens em probabilidades:
  for (s in frq) { frq[s] = frq[s]/tot; }
  printf " %d entradas, total = %d\n", nok, tot > "/dev/stderr";
}

function gera_cadeia(frq,  s,pr)
{
  # Gera uma cadeia com probabilidades especificadas pela tabela {frq}
  pr = rand()*0.9999999;
  for (s in frq)
    { # printf "  %s %.6f %.6f\n", s, pr, frq[s] > "/dev/stderr";
      pr -= frq[s]; 
      if (pr <= 0) { return s; }
    }
  return "***BUG***";
}

function gera_comprimento(tmin,tmed,tmax  )
{
  # Gera um número positivo entre {tmin} e {tmax} com média {tmed}.
  # Por enquanto usa distrinuição uniforme de comprimentos.
  # Deveria usar distribuição exponencial.
 
  # Ajusta {tmin,tmax}:
  if (tmed - tmin < tmax - tmed) { tmax = tmed + (tmed - tmin); }
  if (tmax - tmed < tmed - tmin) { tmim = tmed - (tmax - tmed); }
  
  return tmin + int(rand()*tmax+1-tmin);
}

function arg_error(msg)
{
  printf "** %s\n", msg > "/dev/stderr"; 
  printf "usage: %s\n", usage > "/dev/stderr"; 
  exit 1;
}

function table_error(nome, msg)
{
  printf "%s:%s: ** %s\n", nome, FNR, msg > "/dev/stderr"; 
  exit 1;
}