# -*- cperl -*-

##### General layout of the module #####
#
# Each type of transform that is supported by this module has a plain,
# unthreaded perl entry point the user calls. This entry point makes sure the
# FFTW plan exists (or makes it). Then it calls the THREADED PP function to
# actually compute the transform

# I generate code for up to 10-dimensional FFTs
my $maxrank = 10;

our $VERSION = '0.02.2';

pp_addpm( {At => 'Top'}, scalar `cat README.pod` );

pp_addhdr( '
#include <fftw3.h>
' );


# I want to be able to say $X = fft1($x); rank is required. 'fft()' is ambiguous
# about whether threading is desired or if a large fft is desired. Old PDL::FFTW
# did one thing, matlab does another, so I do not include this function at all




# I define up to rank-10 FFTs. This is annoyingly arbitrary, but hopefully
# should be sufficient
for my $rank (1..$maxrank)
{
  generateDefinitions($rank);
}
pp_export_nothing();

pp_addxs('', `cat compute_plan_template.xs`);
pp_addpm( {At => 'Middle'}, scalar `cat FFTW3_header_include.pm` );


for my $rank (1..$maxrank)
{
  my $shapestr,$rshapestr;
  $shapestr =  sprintf(q{$a->shape->slice('1:%d')->prodover},$rank);
  $rshapestr = sprintf(q{$a->shape->slice('0:%d')->prodover},$rank-1);

  pp_addpm({At => 'Bot'}, <<EOF );
sub fft$rank { __fft_internal( "fft$rank",\@_ ); }
eval( "*PDL::fft$rank = \&fft$rank;" );

sub ifft$rank { my \$a = __fft_internal( "ifft$rank", \@_ ); \$a /= $shapestr; \$a; }
eval( "*PDL::ifft$rank = \&ifft$rank;" );

sub rfft$rank { __fft_internal( "rfft$rank", \@_ ); }
eval( "*PDL::rfft$rank = \&rfft$rank;" );

sub irfft$rank { my \$a = __fft_internal( "irfft$rank", \@_ ); \$a /= $rshapestr; \$a; }
eval( "*PDL::irfft$rank = \&irfft$rank;" );
EOF

  pp_add_exported( "fft$rank", "ifft$rank", "rfft$rank", "irfft$rank" );
}


##########
# Generate the fftn case.  This should probably be done more prettily; for now it's just 
# a springboard that jumps into __fft_internal.
pp_addpm ( {At=> 'Bot'}, <<'EOF' );
sub _rank_springboard {
  my ($name, $source, $rank, @rest) = @_;
  my $inverse = ($name =~ m/^i/);
  my $real    = ($name =~ m/r/);

  unless(defined $rank) {
    die "${name}n: second argument must be the rank of the transform you want";
  }
  $rank = 0+$rank;  # force numeric context
  unless($rank>=1 ) {
    die "${name}n: second argument (rank) must be between 1 and $maxrank";
  }

  my $active_lo = ($real ? 0 : 1);
  my $active_hi = ($real ? $rank-1 : $rank);

  unless($source->ndims > $active_hi) {
    die "${name}n: rank is $rank but input has only ".($active_hi-$active_lo)." active dims!";
  }

  my $out = __fft_internal( $name.$rank, $source, @rest );

  if($inverse) {
    $out /= $out->shape->slice("$active_lo:$active_hi")->prodover;
  }
  return $out;
}

sub fftn    { _rank_springboard( "fft",      @_ ) }
sub ifftn   { _rank_springboard( "ifft",     @_ ) }
sub rfftn   { _rank_springboard( "rfft",  @_ ) }
sub irfftn  { _rank_springboard( "irfft", @_ ) }

EOF
pp_add_exported( map { "${_}fftn" } ('','i','r','ir') );

pp_done();


sub generateDefinitions
{
  my $rank = shift;

  ################################################################################
  ####### first I generate the definitions for the simple complex-complex FFT case
  my $funcname = "__fft$rank";

  # make dimension string 'n0=2,n1,n2,n3,n4...'. The leading 2 is for the
  # (real,imag) complex pair
  my @dims = map {"n$_"} 1..$rank;
  unshift @dims, 'n0=2';
  my $dims_string = join(',', @dims);

  my $code =  `cat template_complex.c`;

  my %pp_def = ( HandleBad    => 0,
                 Pars         => "in($dims_string); [o]out($dims_string);",
                 GenericTypes => [F,D],
                 Code         => $code,
                 OtherPars    => 'IV plan', # comes not from the user, but
                                            # from the pre-fft code

                 # this is a private function so I don't want to create
                 # user-visible documentation or exports
                 Doc          => undef,
                 PMFunc       => ''
               );
  pp_def( $funcname, %pp_def );



  ##################################################################################
  ####### now I generate the definitions for the real-complex and complex-real cases
  my @dims_real    = @dims;
  my @dims_complex = @dims;

  shift @dims_real; # get rid of the (real,imag) dimension for the real numbers
  $dims_complex[1] = 'nhalf'; # first complex dim is real->dim(0)/2+1
  my $dims_real_string    = join(',', @dims_real);
  my $dims_complex_string = join(',', @dims_complex);

  my $code_real          =  `cat template_real.c`;
  $code_real             =~ s/RANK/$rank/ge;
  my $code_real_forward  =  $code_real;
  my $code_real_backward =  $code_real;

  $code_real_forward  =~ s/INVERSE/0/g;
  $code_real_backward =~ s/INVERSE/1/g;

  # forward
  # I have the real dimensions, but not nhalf
  $pp_def{RedoDimsCode} = <<'EOF';
if( $PDL(complex)->ndims <= 1 || $PDL(complex)->dims[1] <= 0 )
  $SIZE(nhalf) = (int)( $PDL(real)->dims[0]/2 ) + 1;
EOF
  $pp_def{Pars} = "real($dims_real_string); [o]complex($dims_complex_string);";
  $pp_def{Code} = $code_real_forward;
  pp_def( "__rfft$rank", %pp_def );

  # backward
  # I have the complex dimensions. Have nhalf, but not n1
  #
  # if we're not given an output, there's an ambiguity. I want
  # int($out->dim(0)/2) + 1 != $in->dim(1),
  # however this could mean that
  #  $out->dim(0) = 2*$in->dim(1) - 2
  # or
  #  $out->dim(0) = 2*$in->dim(1) - 1
  #
  # WITHOUT ANY OTHER INFORMATION, I ASSUME EVEN INPUT SIZES, SO I ASSUME
  #  $out->dim(0) = 2*$in->dim(1) - 2
  $pp_def{RedoDimsCode} = <<'EOF';
if( $PDL(real)->dims[0] <= 0 )
  $SIZE(n1) = 2*$PDL(complex)->dims[1] - 2;
EOF
  $pp_def{Pars} = "complex($dims_complex_string); [o]real($dims_real_string);";
  $pp_def{Code} = $code_real_backward;
  pp_def( "__irfft$rank", %pp_def );
}