#define PROG_NAME "dnafour"
#define PROG_DESC "Fourier analysis of DNA sequences"
#define PROG_VERS "1.1"

/* Copyright © 2003 by the the Fluminense Federal University (UFF).
** See the copyright, authorship, and warranty notice at end of file.
** Last edited on 2007-11-01 20:06:34 by stolfi
*/

#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include <fftw3.h>

int main (int argc, char **argv)
{
  fftw_complex *in, *out;
  fftw_plan plan;

  int nSeqs, nBases;

  /* Leitura de parâmetros de chamada */

  if (argc != 3)
    {
      fprintf(stderr, "Parametros errados!\n");
      fprintf(stderr, "%s N_SEQS TAM_SEQ\n", PROG_NAME);
      return(1);
    }

  nSeqs = atoi(argv[1]);
  nBases = atoi(argv[2]);

  int fMax = nBases/2;  /* Freqüência máxima. */

  /* Aloca de vetor de espectro de potencia e inicializa */
  double pot[fMax+1];
  int i;
  for (i = 0; i <= fMax; i++) { pot[i] = 0; }

  /* Aloca vetores {in,out} para transformada de Fourier: */
  in = fftw_malloc(sizeof(fftw_complex)*nBases);     /* aloca vetor de entrada */
  out = fftw_malloc(sizeof(fftw_complex)*nBases);    /* aloca vetor de saída */
  
  /* Pre-calcula parâmetros para a transformada de Fourier: */
  plan = fftw_plan_dft_1d(nBases, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
  
  /* Fator de normalização para o espectro de potência: */
  double norm = 1.0/nBases;

  /* Calcula FFT para as sequencias e acumula espectros de potência: */
  int s;
  for(s = 1; s <= nSeqs; s++) /* repete para qtdd de sequencias */
    {
      /* Lê um trecho de DNA, converte, guarda no vetor {in}: */
      int nRead = 0;
      while(nRead < nBases)
	{ 
          int base = getc(stdin);
          if (base != '\n')
            {
              /* Mapeia bases {A,T,C,G} para números complexos {±1,±I}: */
              switch(base)
                { 
                case 'A': in[nRead][0] = +1; in[nRead][1] = 00; break;
                case 'T': in[nRead][0] = -1; in[nRead][1] = 00; break;
                case 'C': in[nRead][0] = 00; in[nRead][1] = +1; break;
                case 'G': in[nRead][0] = 00; in[nRead][1] = -1; break;
                default:
                  fprintf(stderr, "Ivalid base (%c)\n", base);
                  return 1;
                }
              nRead++;
            }
	}


      /* Calcula transformada de Fourier */
      fftw_execute(plan);

      /* Calcula espectro de potência */
      /* Stolfi: "Lembre-se que no espectro de potência deve-se somar os
         módulos ao quadrado de coeficientes de mesma freqüência, isto é
         {|out[i]|^2 + |out[N-i]|^2} , para {i = 1..N/2-1}; sendo que 
         {out[0]} (freqüência 0) e {out[N/2]} (freqüência {N/2}) 
         ficam sem parceiro. Portanto o espectro de potência 
         tem elementos {pot[0..N/2]}, inclusive ambas as pontas."
      */
      
      double poti, potj; /* {= |out[i]|^2, |out[j]|^2} */
      
      /* Freqüência 0: */
      i = 0;
      poti = out[i][0]*out[i][0] + out[i][1]*out[i][1];
      pot[i] += norm * poti;
      
      /* Freqüência {1..(N/2)-1}: */
      for (i = 1; i <= fMax-1; i++)
	{
          int j = nBases-i;
          poti = out[i][0]*out[i][0] + out[i][1]*out[i][1];
          potj = out[j][0]*out[j][0] + out[j][1]*out[j][1];
          pot[i] += norm * (poti + potj);
	}
        
      /* Freqüência {N/2}: */
      i = fMax;
      poti = out[i][0]*out[i][0] + out[i][1]*out[i][1];
      pot[i] += norm * poti;

    }

  /* Calcula média dos espectros de potência: */
  for (i = 0; i <= fMax; i++)
    {
      fprintf(stdout, "%6d %12.4f\n", i, pot[i]/nSeqs);
    }

  fclose(stdout);

  fftw_destroy_plan(plan);
  fftw_free(in); fftw_free(out);
  return 0;
}

/* COPYRIGHT, AUTHORSHIP, AND WARRANTY NOTICE:
**
**   Copyright © 2003 by the the Fluminense Federal University (UFF).
** 
** Created 2003 by Luciana Pessôa, UFF.
** Heavily modified 2004-2005 by Jorge Stolfi, Unicamp.
**
** Permission to use, copy, modify, and redistribute this software and
** its documentation for any purpose and without fee is hereby
** granted, provided that: (1) the copyright notice at the top of this
** file and this copyright, authorship, and warranty notice is retained
** in all derived source files and documentation; (2) no executable
** code derived from this file is published or distributed without the
** corresponding source code; and (3) these same rights are granted to
** any recipient of such code, under the same conditions.
** This software is provided "as is", WITHOUT ANY EXPLICIT OR IMPLICIT
** WARRANTIES, not even the implied warranties of merchantibility and
** fitness for a particular purpose. END OF NOTICE.
*/