use strict;
use Carp;

pp_bless('PDL::Fit::Levmar::Func');  # make the functions generated go into our namespace     
pp_addpm({At=>'Top'},<<'======EOD_Documentation======');
=head1 NAME

PDL::Fit::Levmar::Func - Create model functions for Levenberg-Marquardt fit routines

=head1 DESCRIPTION

This module creates and manages functions for use with the
Levmar fitting module.  The functions are created with a
very simple description language (that is mostly C), or are
pure C, or are perl code. For many applications, the present
document is not necessary because levmar uses the
Levmar::Func module transparantly.  Therefore, before
reading the following, refer to L<PDL::Fit::Levmar>.

=head1 SYNOPSIS

    use PDL::Fit::Levmar::Func;

   $func = levmar_func(  FUNC =>  # define and link function

 '
  function gaussian
  loop
  x[i] = p[0] * exp(-(t[i] - p[1])*(t[i] - p[1])*p[2]);
  end function

  jacobian jacgaussian
  double arg, expf;
  loop
  arg = t[i] - p[1];
  expf = exp(-arg*arg*p[2]);
  d0 = expf;
  d1 = p[0]*2*arg*p[2]*expf;
  d2 = p[0]*(-arg*arg)*expf;
  end jacobian 
  ' 
  );

  $hout = PDL::Fit::Levmar::levmar($p,$x,$t,  # fit the data
               FUNCTION => $func,
	 );

=head1 FUNCTIONS

=head2 levmar_func()

=for ref

This function creates and links a function, ie, takes a function definition
and returns a function (object instance) ready for use by levmar.
see PDL::Fit::Levmar::Func for more information.

=for usage

 levmar_func

=for example


=for options

=head2 call()

=for ref

Call (evaluate) the fit function in a Levmar::Func object.

=for usage

 use PDL::Fit::Levmar::Func;

 $Gh = levmar_func(FUNC=>$fcode);

 $x = $Gh->call($p,$t);

 Here $fcode is a function definition (say lpp code). (This does
 not currently work with perl subroutines, I think. Of course,
 you can call the subroutine directly. But it would be good
 to support it for consistency.)
 
 $p is a pdl or ref to array of parameters.
 $t is a pdl of co-ordinate values.
 $x is the fit function evaluated at $p and $t. $x has the
 same shape as $t.

=for example

 $Gf = levmar_func( FUNC => '
       function	   
       x = p0 * t * t;
     ');

 print $Gf->call([2],sequence(10))  , "\n";
 [0 2 8 18 32 50 72 98 128 162]

=cut 


use strict;
use Carp;
use Config;
use File::Path;
use File::Spec::Functions;

sub munlink;

#use PDL::NiceSlice;
use PDL::Core ':Internal'; # For topdl()

use vars ('$MKOBJ', '$MKSO', '$OBJ_EXT' , '$SO',
	  '$LEVFUNC', '$JLEVFUNC', '$LPPEXT' );

# These are pointers to the C functions that wrap the user's
# perl fit functions.
#$LEVFUNC = get_perl_func_wrap();
#$JLEVFUNC = get_perl_jac_wrap();

$OBJ_EXT = $Config{obj_ext};
$SO = $Config{so};
$LPPEXT = '.lpp';
# string  to compile .c to .o 
$MKOBJ = $Config{cc} . " -c " . $Config{optimize} . " " .
    $Config{cccdlflags} . " ";

# string  to compile .o to .so 
$MKSO = $Config{ld} . " " . $Config{lddlflags} ." %o -o %s ";


======EOD_Documentation======

pp_addpm({At=>'Bot'},<<'======EOD_Authors======');

=head1 AUTHORS

Copyright (C) 2006 John Lapeyre.
All rights reserved. There is no warranty. You are allowed
to redistribute this software / documentation under certain
conditions. For details, see the file COPYING in the PDL
distribution. If this file is separated from the PDL distribution,
the copyright notice should be included in the file.

=cut

======EOD_Authors======

pp_addhdr('
#include <stdio.h>
#include <stdlib.h>
#include "func.h"

');

pp_add_exported('', 'levmar_func');
pp_addpm(<<'======EOD_new======');

=head2 new

=for usage

  my $fitfunc = new PDL::Fit::Levmar::Func(...);
  or
  my $fitfunc = levmar_func(...);

=for ref

Generic constructor 

=cut

$PDL::Fit::Levmar::Func::SRC_EXT = ".c";
$PDL::Fit::Levmar::Func::SO_EXT = ".so";
$PDL::Fit::Levmar::Func::OBJ_EXT = ".o";

sub deb { print STDERR $_[0], "\n";}

##########################################
# Methods follow

sub new {
  my($class) = shift;
  my($opts) = shift;
  if(ref $opts ne 'HASH') { # change opts to hash of opts
    $opts = defined $opts ? {$opts,@_} : {} ;
  }
  my $self = { # following are defaults
      MKOBJ => $MKOBJ,
      MKSO => $MKSO,
      NOCLEAN => undef,
      FUNC => undef,
      JFUNC => undef,
#      TYPE => 'double',
      NOSO => undef,
      CSRC => undef,
      DIR =>  "./tempcode" , # where the .so is built
      TESTSYNTAX => undef,
      LIBHANDLE => undef,
      
  };
  foreach ( keys %$opts ) { # replace defaults
      if ( exists $self->{$_} ) {
	  $self->{$_} = $opts->{$_};
      }	  
      else { die "Levmar::Func::new: Unrecognized option to Levmar::Func '$_'";}
  }
  return bless $self,$class;
}

# all these attempts at dying and croaking fail.
# program does not quit and error messages are not printed
# dont know why. It must be some bad bug in the xs code.
# If I just print the error message, its fine.
# I think it has to do with DESTROY. return does not work either
sub DESTROY {
    my $self = shift;
    my $cnt = 0;
    return if defined $self->{NOSO};
    if ( defined $self->{LIBHANDLE}) { 
	my ($ret, $err) = close_shared_object_file($self->{LIBHANDLE});
	if ( $err ne '' or $ret != 0) {
	    warn "Levmar::Func::DESTROY: Error closing function library:$ret '$err'";
	}
	warn "Error closing function library:$ret '$err'" if  $ret != 0;
	if ( defined $self->{JPOINTER} ) { # call close twice because we got two symbols
	    ($ret, $err) = close_shared_object_file($self->{LIBHANDLE});
	    warn "Error closing function library:$ret '$err'" if $ret !=0 ;
	    if ( $err ne '' or $ret != 0 ) {
		warn "Levmar::Func::DESTROY: Error closing function library:$ret '$err'";
	    }
	}
	if ( defined $self->{SJPOINTER} ) { # call close twice because we got two symbols
	    ($ret, $err) = close_shared_object_file($self->{LIBHANDLE});
	    warn "Error closing function library:$ret '$err'" if $ret !=0 ;
	    if ( $err ne '' or $ret != 0 ) {
		warn "Levmar::Func::DESTROY: Error closing function library:$ret '$err'";
	    }
	}
	if ( defined $self->{SFPOINTER} ) { # call close twice because we got two symbols
	    ($ret, $err) = close_shared_object_file($self->{LIBHANDLE});
	    warn "Error closing function library:$ret '$err'" if $ret !=0 ;
	    if ( $err ne '' or $ret != 0 ) {
		warn "Levmar::Func::DESTROY: Error closing function library:$ret '$err'";
	    }
	}
    }
    if (not defined $self->{NOCLEAN} and defined $self->{SONAME} ) {
	foreach my $k ( qw ( SONAME  )) { # clean so file
	    my $fn = $self->{$k};
	    if (-e $fn) {
		if ( not -f $fn ) {
		    warn "Levmar::Func::DESTROY: Refusing to remove non-file '$fn'";
		}
		else {
		    $cnt  = munlink $fn;
		    if ( $cnt < 1 ) {
			warn "Levmar::Func::DESTROY: failed removing '$fn' (lib probably open)";
		    }
		    else {
#			warn "Levmar::Func::DESTROY: SUCCESS removing '$fn'";
		    }
		}
	    }
	}
    }
}

sub levmar_func { my $self = new PDL::Fit::Levmar::Func(@_);
		  my @err;
		  if ( defined $self->{FUNC} and $self->{FUNC} =~ /CODE/ ) {
		      $self->{NOSO} = 1; # bad way to handle this
		      $self->{FPOINTER} = 0; # test for this in levmar
		  }
		  else {
		      if ( not defined $self->{FUNC}  and not defined $self->{CSRC} ) {
			  die "Levmar::Func::lemar_func: neither FUNC nor CSRC defined";
		      }
		      @err = $self->make_ccode();
		      return @err if $err[0] == -1;
		      $self->make_build_names();
		      $self->write_ccode();
		      if ( not $self->make_shared_object_file() ) {
			  die "Levmar::Func::levmar_func: Compilation failed.";
		      }
		      if ( not $self->load_fit_library() ) {
			  die "Levmar::Func::levmar_func: Loading dynamic library failed.";
		      }
		      $self->clean_files() unless defined $self->{NOCLEAN};
		  }
		  if ( defined $self->{JFUNC}  ) { # probably won't get called
		      if ( $self->{JFUNC} =~ /CODE/ ) {
			  $self->{NOSO} = 1;
		      }
		      else {
			  die "Levmar::Func::levmar_func option JFUNC must be a ref to a perl sub";
		      }
		  }
		  return $self;
}


sub make_build_names {
    my $self = shift;
    die "Levmar::Func::levmar_func:  no name for function"
	unless defined $self->{NAME};
    my $q = $self->{DIR};
    if ( not -d $q ) {
	die "Levmar::Func::levmar_func: Can't make directory (folder) '$q'"
	    unless File::Path::mkpath( $q );
    }
    $self->{FNAME} = $self->{NAME} unless exists $self->{FNAME}; # base file name
    my $n = $self->{FNAME};
    my $srcn = catfile( $q , $n . $PDL::Fit::Levmar::Func::SRC_EXT);
    my $son = catfile($q ,$n . "." .  $SO);
    my $objn = catfile($q , $n . $OBJ_EXT);
    $self->{SRCNAME} = $srcn;
    $self->{SONAME} = $son;
    $self->{OBJNAME} = $objn;
}



# change this to identity, can remove sometime
sub mfqp {
    my ($self,$key) = @_;
    return  $self->{$key};
}

# Write c source file from function definition.
# Entire definition is given as a string, or the
# filename , ending in $LPPEXT is given
# Copying the source strings around, inefficient. should use
# refs
sub make_ccode {
    my $self = shift;
    my @err;

    if ( not defined $self->{CSRC} ) { # is FUNC really c code?
	if (not $self->{FUNC} =~ /$LPPEXT$/o        # is not a $LPPEXT file
                and ( $self->{FUNC} =~ /\.c$/ or # is a .c file
                                                 # doesnt have function def
 	 not $self->{FUNC} =~ /(^|\n)\s*function(\s+\w+|\s*)\n/) ) {
	    $self->{CSRC} = $self->{FUNC};
	}
    }
    if ( defined $self->{CSRC} ) {
	if ( $self->{CSRC} =~ /\.c$/ ) {
	    my $cf = $self->{CSRC};
	    local (*CHAND);
	    local $/; 
	    open CHAND, "<$cf" or die "Levmar::Func::make_code : Can't open C src file '$cf'";
	    $self->{CCODE} = <CHAND>;
	    close(CHAND);
	}
	else {  # inefficient, but it works
	    $self->{CCODE} = $self->{CSRC}; 
	}
	my ($funcname, $jacname) = $self->get_funcs_from_c();
	my $floatcode = $self->{CCODE};
#	if ( /^\s*void\s+([\w]+)/ )
#	deb $floatcode;
	$floatcode =~ s/$funcname/s$funcname/;
	$floatcode =~ s/$jacname/s$jacname/ if defined $jacname and $jacname ne ''; # gjl Dec 2011
#	deb $floatcode;
	$self->{CCODE} = 
"
#define FLOAT double
" .  $self->{CCODE} .
"
#undef FLOAT
#define FLOAT float
" . $floatcode;
	$self->{NAME} = $funcname unless exists $self->{NAME};
	$self->{JACNAME} = $jacname unless exists $self->{JACNAME};
	return (1); # don't use FUNC even if its there.
    }
    elsif ( defined $self->{FUNC} ) {
	my $st = $self->{FUNC};
	if ($st =~ /$LPPEXT$/o) {  # file, not string
	    local (*FUNCHAND);
	    local $/; 
	    open FUNCHAND, "<$st" or die "Levmar::Func::make_ccode : Can't open def file '$st'";
	    $st = <FUNCHAND>;
	    close(FUNCHAND);
	}
	my ($funcname,$jacname,$ccode);
	my @ret = generate_C_source($self, $st);
	if ( $ret[0] eq '-1' ) { # stop warnings with ==
	    my $r;
#	    ($r, $funcname, $jacname, $code, @err) = @ret; 
	    return @ret;
	}
	else {
	    ($funcname,$jacname,$ccode) = @ret;
	}
	$self->{CCODE} = $ccode;
	$self->{NAME} = $funcname unless exists $self->{NAME};
	$self->{JACNAME} = $jacname unless exists $self->{JACNAME};
	return (1);
    }
    return (-1, "Can't make C code with no defintion or source file");
}

# jacobian function name **must** start with 'jac' if using pure c code
sub get_funcs_from_c {
    my $self = shift;
    my @lines = split ("\n",$self->{CCODE});
    my ($fname, $jacname) = (undef, undef);
    foreach (@lines ) {
	if ( /^\s*void\s+([\w]+)/ ) { # prototype begins with 'void' return type
	    my $fn = $1;
	    if ($fn =~ /^jac/ ) {
		$jacname = $fn;
	    }
	    else {
		$fname = $fn;
	    }
	}
    }
    if ( not defined $fname or $fname eq '' ) { # ok if jacname == undef
	die "Levmar::Func::get_funcs_from_c : Can't find function name from c code";
    }
    return ($fname, $jacname);
}

sub write_ccode {
    my $self = shift;
    local (*CCODEH);
    return unless exists $self->{CCODE};
    if ( not defined $self->{NAME} ) {
	die "Levmar::Func::write_ccode : Can't write C code for function with no name";
    }
    my $srcn = $self->{SRCNAME};
    open CCODEH, ">$srcn" or die "Levmar::Func::write_ccode : Can't open C source file '$srcn' for writing";
    print CCODEH $self->{CCODE};
    close CCODEH;
}

sub make_shared_object_file {
    my $self = shift;
    my $ret;
    my $srcn1 = $self->{SRCNAME};
    my $fname1 = $self->{FNAME};
    my $i = 0;
    my $max = 20; # allow $max renamed files
    while ( -f $self->{SONAME} ) { # this is important to avoid user side bugs
	$i++;
	my $so = $self->{SONAME};
	$self->{FNAME} = $fname1 . $i;
	$self->make_build_names();
#	warn "Levmar::Func::make_shared_object_file: $so exists. Trying " . $self->{SONAME};
	last if $i > $max;
    }
    my $srcn = $self->{SRCNAME};
    my $objn = $self->{OBJNAME};
    my $sobjn = $self->{SONAME};
    my $fname = $self->{FNAME};
    if ( $i > $max ){
	die 'Levmar::Func::make_shared_object_file: Too many renamed files in '. $self->{DIR} . ' try cleaning.';
    }
    elsif ($i > 0) {
	if ( rename($srcn1,$srcn) == 0 ) {
	    die "Levmar::Func::make_shared_object_file: Faile to rename $srcn1 to $srcn";
	}
    }
    my $mkobjstr = $self->{MKOBJ} ." -o $objn  $srcn";
    $ret = do_system( $mkobjstr); # compile c to o
#    deb "Levmar::Func::make_shared_object_file: compile returned '$ret'";
    if ( $ret != 0) { # zero is success
	warn "Levmar::Func::make_shared_object_file: Compiling '$srcn' failed.";
	return 0;
    }
    $self->{MKOBJSTR} = $mkobjstr;
    my $sostr = $self->{MKSO};
    $sostr =~ s/\%s/$sobjn/g;
    $sostr =~ s/\%o/$objn/g;
    $ret = do_system( $sostr ); #  make so from o
#    deb "Levmar::Func::make_shared_object_file: '$sostr'  make so from o returned '$ret'";   
    if ( 0 != $ret ) {
	warn "Levmar::Func::make_shared_object_file: Making shared library '$sobjn' failed.";
	return 0;
    }
    $self->{SOSTR} = $sostr;
}

#### Some methods to return information
# get compiler strings
sub get_cc_str {
    my $self = shift;
    return ($self->{MKOBJSTR}, $self->{SOSTR});
}

sub get_fpoint {
    my $self = shift;
    $self->{FPOINTER};
}

sub get_jacpoint {
    my $self = shift;
    $self->{JPOINTER};
}

sub get_fit_pointers {
    my $self = shift;
    ($self->{FPOINTER},
     $self->{SFPOINTER},
     $self->{JPOINTER},
     $self->{SJPOINTER});
}

sub get_file_names {
    my $self = shift;
    return ($self->mfqp('SRCNAME'),
	    $self->mfqp('OBJNAME'),
	    $self->mfqp('SONAME'));
}

# remove generated  o, and C files
# Well, don't remove so file.
sub clean_files {
    my $self = shift;
    my $cnt = 0;
    my $srcn = $self->mfqp('SRCNAME'); # fully qualified names
    my $objn = $self->mfqp('OBJNAME');
#    foreach my $k ( qw ( OBJNAME SONAME )) {
    foreach my $k ( qw ( OBJNAME  )) { # soname needs to be there
	my $fn = $self->mfqp($k);
	if (-e $fn) {
	    if ( not -f $fn ) {
		warn "Levmar::Func::clean_files : Object file '$fn' not a file. Not deleting";
	    }
	    else {
		$cnt  = munlink $fn;
		if ( $cnt < 1) {
		    warn "Levmar::Func::clean_files : cnt $cnt: Failed to delete object file '$fn'";
		}
	    }
	}
    }
    if ( -e $srcn ) {
	if ( not -f $srcn ) {
	    warn "Levmar::Func::clean_files: Source file '$srcn' not a file. Not deleting";
	}
	local (*SRCH);
	open SRCH, "<$srcn" or die "Levmar::Func::clean_files : Can't open C source file '$srcn' for reading";
	while (<SRCH>) {
	    next if /^\s*$/;
	    if (/\/\* This file Automatically/) {
		$cnt = 0;
		$cnt = munlink $srcn;
		if ( $cnt < 1) {
		    warn "Levmar::Func::clean_files: cnt $cnt: Failed to delete source file '$srcn'";
		}
	    }
	    last; # only check up to first non-blank line
	}
	close(SRCH);
    }
}

# Improve effiency after this is well tested by getting all four
# symbols at once.

sub load_fit_library {
    my $self = shift;
    my $so = $self->mfqp('SONAME');
    my $fn = $self->{NAME};
    my $jn = $self->{JACNAME};
    die "Levmar::Func::load_fit_library : Can't find dynamic library file '$so'" unless -e $so;
    my ($func_pointer, $jac_pointer, $lib_handle, $error_message)=
	(0,0,0,0);
    my ($sfunc_pointer, $sjac_pointer) =    # single precision
	(0,0);
    ($func_pointer, $lib_handle, $error_message) =
	open_and_find($so,$fn);
    ($sfunc_pointer, $lib_handle, $error_message) =  # single precision
	open_and_find($so,'s'.$fn);
    if ( $error_message ne ''  ) { 
	warn "Levmar::Func::load_fit_library: Failed to load dynamic library '$so'\n".
          " and find fit func '$fn'.\n dlopen:'$error_message'";
	return 0;
    }
    if ( defined $jn and $jn ne '') {
	($jac_pointer, $lib_handle, $error_message) =
	    open_and_find($so,$jn);
	if ( $error_message ne ''  ) { 
	    warn "Levmar::Func::load_fit_library: Failed to load dynamic library '$so' and find jacobian '$jn'.\n".
		" dlopen:'$error_message' ";
	}
	else {
	    $self->{JPOINTER} = $jac_pointer;
	}
	($sjac_pointer, $lib_handle, $error_message) =
	    open_and_find($so, 's' . $jn);
	if ( $error_message ne ''  ) { 
	    warn "Levmar::Func::load_fit_library: Failed to load dynamic library '$so' and find jacobian '$jn'.\n".
		" dlopen:'$error_message' ";
	}
	else {
	    $self->{SJPOINTER} = $sjac_pointer;
	}

    }
    $self->{FPOINTER} = $func_pointer;
    $self->{SFPOINTER} = $sfunc_pointer;
    $self->{LIBHANDLE} = $lib_handle;
}


# End of methods
# (Actually there are some more methods near the bottom
#  that are close to the corresponding pp_defs)
########################################
# Non-methods below here

# should return 0 on success
sub do_system {
    my $c = shift;
#    print STDERR $c,"\n";
    system $c;
}

sub munlink {
    my $f = shift;
    my $cnt;
    my $sf = "'$f'";
    if ( not -e $f ) {
	deb "munlink: $sf does not exist";
	return 0;
    }
    if ( not -f $f ) {
	deb "munlink: $sf is not a file. Not removing.";
	return 0;
    }
    $cnt = unlink $f;
    if (not $cnt > 0) {
	deb "munlink: Can't unlink exisiting file $sf";
    }
    return $cnt;
}

# Turn a function description into C source
# !care! I disabled NiceSlice, because it erroneously saw
# slicing in this piddle-less routine. But I was not using it
# anyway. (parser is rewritten, maybe not a problem anymore.)
sub generate_C_source {
    my ($self, $fdef) = @_;
    my $funcname = ''; # this one is returned
    my $jacname = ''; # this one is returned
    # Do this first now send t -> t[i], x -> x[i], p3 -> p[3]
    # Probably a better way than doing this in a loop
    # but overlapping matchin patterns do not work with global flag
    # t --> t[i]

    # I tried to substitutions on $fdef vs on each line. Not much difference.
    my $defs = {};
    my ($lines,$def, $loopf, $noloopf, $name);
    my (@flines) = split(/\n/, $fdef);
    $lines = [];
    my $prefunc = $lines;# lines before any func def
    foreach my $oneline (@flines) { #break code into functions
	next if $oneline =~ /^\s*end\s+(function|jacobian)\s*/;
	if ( $oneline =~ /^\s*(function|jacobian)\s*(\w*)\s*$/ ) { # func declaration
	    $name = $2;
	    $name = "func" unless $name ne '';
	    my $type = $1;
	    $name = "jac" . $name if $type eq 'jacobian' and not $name =~ /^jac/ ;
	    $def = $defs->{$name} = {};
	    $def->{type} = $type;
	    $lines = [];
	    $def->{loop} = $lines;
	    $noloopf = 0;
	    $loopf  = 0;
	    next;
	}
	if ( $oneline =~  /^\s*noloop/ ) {  # don't want any implicit loop
	    $def->{preloop} = $lines;
	    $def->{loop} = undef;
	    $noloopf = 1;
	    next;
	}
	elsif ( $oneline =~  /^\s*loop/ ) {  # start implicit loop
	    $def->{preloop} = $lines;
	    $lines = [];
	    $def->{loop} = $lines;
	    $loopf = 1;
	    next;
	}
	else {
	    push @$lines, "" .$oneline ."\n";
	}
    }

    
    if ( defined $self->{TESTSYNTAX} ) {
	my $st = '';
        if (defined @$prefunc > 0) {
	    $st .= ":PREFUNC: lines:" . scalar(@{$prefunc}). "\n";
	    	 $st .= join("\n",@{$prefunc}). "\n";
	}	    
	foreach my $name ( keys %$defs ) {
	    my $d = $defs->{$name};
	    $st .=  "NAME: $name, TYPE: " . $d->{type};
	    if (defined $d->{preloop}) {
		$st .=  ":PRELOOP: lines:" . scalar(@{$d->{preloop}}). "\n";
		$st .=  join("\n",@{$d->{preloop}}). "\n";
	    } 
	    if (defined $d->{loop}) {
		$st .=  ":LOOP: lines:" . scalar(@{$d->{loop}}). "\n";
		$st .=  join("\n",@{$d->{loop}}). "\n";
	    }
	}
	$self->{SYNTAXRESULTS} = $st;
    }

    my $end_function = "\}\n\n";
    my $end_loop = "  \}\n";
    my $start_loop = { 'jacobian' =>  "  for(i=j=0; i<n; ++i){\n",
		       'function' =>  "  for(i=0; i<n; ++i){\n" };

#    my $dtype = $self->{TYPE};
    my $dtype = 'FLOAT';
    my $start_function = {
	'jacobian' =>  "void JACFUNNAME($dtype *p, $dtype *d, int m, int n, $dtype *t)\n\{\n".
	    "  int i,j;\n",
#	    "  $dtype *t = ($dtype *) data;\n",
       'function' =>  "void FITFUNNAME($dtype *p, $dtype *x, int m, int n, $dtype *t)\n\{\n".
                 "  int i;\n" };
#             "  $dtype *t = ($dtype *) data;\n" };
    my $kode = 
	"/* This file Automatically created and deleted. Do not edit. */\n".
	"#include<math.h>\n#include <stdio.h>\n\n" .
	join('',@$prefunc);


    my $kode1 = '';

    foreach my $name ( keys %$defs ) {
	my $d = $defs->{$name};
	my ($loop, $preloop, $loopcode);
	my $type = $d->{type};
	$funcname = $name if $type eq 'function';
	$jacname = $name if $type eq 'jacobian';
	if ( defined $d->{loop} ) {
	    $loopcode = join("", @{$d->{loop} });
	    if ( $loopcode =~ /^\s*$/ ) { # loop code is whitespace
		$loopcode = undef;    # so don't print it
	    }
	    else {
		$loop = $d->{loop};
	    }
	}
	$preloop = $d->{preloop} if defined $d->{preloop};

#	$kode1 .= sprintf($start_function->{$type},$name);
	$kode1 .= sprintf($start_function->{$type}, '');
	if (defined $preloop ) {
	    foreach (@$preloop) {
		# p5 --> p[5]
		while ($_ =~ s/([^\w]|^)p(\d+)([^\w\[])/$1p[$2]$3/g){}
		if ( $type eq 'function' ) {
		    # x7 -- > x[7]
		    while($_ =~ s/([^\w]|^)x(\d+)([^\w\[])/$1x[$2]$3/g){}
		}
		elsif ( $type eq 'jacobian' ) {
		    # d3[6] -> d[ 3+ (n-1)*6], no! [3+m*6] is correct
		    while($_ =~ /([^\w]|^)d(\d+)\[(\d+)\]/ ) {
			my $i = "$2+m*$3" ;
			$_ =~ s/([^\w]|^)d\d+\[\d+\]/$1d[$i]/;
		    }
		}
	    }
	    $kode1 .= join('',@$preloop);
	}
	if ( defined $loop ) {
	    $kode1 .= $start_loop->{$type};
	    foreach ( @$loop ) {
		# t -- > t[i]
		while($_ =~ s/([^\w]|^)t([^\w\[])/$1t[i]$2/g){}
		# p5 --> p[5]
		while ($_ =~ s/([^\w]|^)p(\d+)([^\w\[])/$1p[$2]$3/g){}
		if ( $type eq 'function' ) {
		    # x -> x[i]
		    while($_ =~ s/([^\w]|^)x([^\w\[])/$1x[i]$2/g){}
		}
		elsif ( $type eq 'jacobian' ) {
		    if ( /^\s*d\d+\[?i?\]?(.+)/ ) {
			my $line = $1;
			$line  =~ s/\s*//; # strip space
			$_ =   "   d[j++] $line\n";
		    }
		}
		else {
		    warn "Levmar::Func::generate_C_source: Function type neither 'function' nor 'jacobian'";
		}
	    }
	    $kode1 .= join('',@$loop) . $end_loop;
	}
	$kode1 .= $end_function;
    } #foreach my $name ..
    $kode .= "
#define FLOAT double
#define JACFUNNAME $jacname
#define FITFUNNAME $funcname
" .  $kode1 .
"
#undef FLOAT
#undef JACFUNNAME
#undef FITFUNNAME
#define FLOAT float
#define JACFUNNAME s$jacname
#define FITFUNNAME s$funcname
" . $kode1;
    return ($funcname,$jacname,$kode);
}  # end sub generate_C_source {

# load shared library with fit function and find a symbol
# (pointer to fit function or jacobian in our case.)
# $lib_name : name of library file , eg "./func.so"
# $func_name : name of function, eg "square"
# $func_pointer : value of C level pointer to $func_name, cast to int
# $lib_handle : same, but for pointer to library
sub open_and_find {
    my ($lib_name, $func_name) = @_;
    my $lib_handle = 0;
    my $func_pointer = 0;
    my $nchar = 1000; # avoid buffer overun by sending this number to memccpy
    my $error_message = "." x $nchar;
    call_open_and_find($lib_name, $func_name, $lib_handle,
		       $func_pointer,  $error_message , $nchar);
    return ($func_pointer, $lib_handle, $error_message);
}

sub close_shared_object_file {
    my ($lib_handle) = @_;
    my $nchar = 1000; # avoid buffer overun by sending this number to memccpy
    my $error_message = "." x $nchar;
    my $ret = 0;
    call_close_shared_object_file( $lib_handle, $ret, $error_message, $nchar);
    return ($ret, $error_message);
}

#-----------------------------------------
# Some methods to call the functions

# evaluate the function.
# Given parameters p and coordiante pdl t,
# return x(t;p) with x->dims == t->dims.
# The functions don't really have to have
# this form, but this is what we assume for
# this call.
sub call {
    my $self = shift;
    my ($p,$t) = @_; # pdls
    my $x = zeroes($t);
    if ($self->{FPOINTER} == 0) {
	warn "Levmar::Func::call : No function loaded.";
	return PDL->null;
    }
    _callf(topdl($p),$t,$x, $self->{FPOINTER});
    return($x);
}

# call jacobian and get back 2-d pdl
# return jac(m,n) = jacfunc(p(m),t(n))
sub jac_of_t {
    my $self = shift;
    my ($p,$t) = @_; # pdls
    if ($self->{JPOINTER} == 0) {
	warn "Levmar::Func::jac_of_t: No jacobian function loaded.";
	return PDL->null;
    }
    my @dims = dims($t);
    my $jac = zeroes($p->type, $p->nelem,@dims);
    return 0 unless ref($jac) =~ /PDL/;
    _callj($p,$t,$jac, $self->{JPOINTER});
    return($jac);
}

# treat array of derivatives as 1-d, just at the
# supplied c routines do. For debugging.
# return jac(m*n) = jacfunc(p(m),t(n))
sub jac_of_t1 {
    my $self = shift;
    my ($p,$t) = @_; # pdls
    if ($self->{JPOINTER} == 0) {
	warn "Levmar::Func::jac_of_t1: No jacobian function loaded.";
	return PDL->null;
    }
    my @dims = dims($t);
    foreach (@dims) { $_ *= $p->nelem}
    my $jac = zeroes($p->type, @dims);
    return 0 unless ref($jac) =~ /PDL/;
    _callj1($p,$t,$jac, $self->{JPOINTER});
    return($jac);
}

======EOD_new======

pp_def ( '_callf' ,
       Pars => ' double p(m); double t(nt); double [o]x(n);',
       OtherPars => ' IV func; ',
       GenericTypes => ['D'],  Doc => undef,
       Code => ' void (*fp)(double *, double *, int, int, void *);
                 fp = (void (*)(double*, double*, int, int, void*)) $COMP(func);
                 fp($P(p),$P(x),$SIZE(m),$SIZE(n),(void *)$P(t));
               ');

# This treats d as 2d pdl
pp_def ( '_callj' ,
       Pars => ' double p(m); double t(nt); double [o]x(m,n);',
       OtherPars => ' IV func; ',
       GenericTypes => ['D'], Doc => undef,
       Code => ' void (*fp)(double *, double *, int, int, void *);
                 fp = (void (*)(double*, double*, int, int, void*)) $COMP(func);
                 fp($P(p),$P(x),$SIZE(m),$SIZE(n),(void *)$P(t));
               ');

# d is longer 1d pdl
pp_def ( '_callj1' ,
       Pars => ' double p(m); double t(nt); double [o]jac(jm);',
       OtherPars => ' IV func; ',
       GenericTypes => ['D'], Doc => undef,
       Code => ' void (*fp)(double *, double *, int, int, void *);
                 fp = (void (*)(double*, double*, int, int, void*)) $COMP(func);
                 fp($P(p),$P(jac),$SIZE(m),$SIZE(nt),(void *)$P(t));
               ');

# Following routines are in getsymbol.c and do not invovle piddles
# note that pointers to ints are passed for library and function pointers
# The routines called use them as pointers to pointers to void (ie for library
# handle and functions).
# I don't know how portable this is. It relies on ints and pointers being the
# same size.

pp_addxs ( '', '

MODULE = PDL::Fit::Levmar::Func PACKAGE = PDL::Fit::Levmar::Func


void
call_open_and_find( lib_name, func_name, lib_handle, func_pointer,  error_message, nchar )
  char * lib_name
  char * func_name
  void *  lib_handle
  void *  func_pointer
  char * error_message
  int  nchar
  CODE:
  _open_and_find( lib_name, func_name, &lib_handle, &func_pointer,  error_message, nchar );
  OUTPUT:
  lib_handle
  func_pointer
  error_message

void
call_close_shared_object_file ( lib_handle, rval , error_message, nchar )
   void *  lib_handle
   int  rval
   char * error_message
   int   nchar
   CODE:
   _close_shared_object_file ( &lib_handle, &rval , error_message, nchar);
   OUTPUT:
   rval
   error_message


');