/* This file is to be included in Levmar.xs
 *
 * John Lapeyre
 */

#include "EXTERN.h"
#include "perl.h"
#include "XSUB.h"
#include "pdl.h"
#include "pdlcore.h"
#include <stdio.h>

static Core* PDL; /* Structure holds core C functions */
static SV* CoreSV;       /* Gets pointer to perl var holding core structure */

// following line useless, i guess
// typedef void (*levmarfunc)( double *, double *, int, int, void * ) ;

typedef void (*DelMagic)(pdl *, int param);
static void delete_levmar_pdls(pdl* p, int param);

static void default_magic(pdl *p, int pa) { p->data = 0; }
static pdl* pdl_wrap(void *data, int datatype, PDL_Long dims[],
		     int ndims, DelMagic delete_magic, int delparam);
static SV *getref_pdl(pdl *it);

/* We use a struct holding data for the wrapper functions to
 * perl fit code, LEVFUNC, and JLEVFUNC. liblevmar uses fit
 * functions whose the last argument is a pointer to
 * unpecified data. When fitting data it is normally used
 * for the ordinate data (t). We pass a DFP struct instead.
 * These functions are all re-entrant, i hope. The DFP
 * struct holds pdls which are created once at the start of
 * a fit and destroyed when it is finished. Data is not put
 * in the pdls here. The fit routines send different arrays
 * on different calls, so the pdl pointer to the data must
 * be set in the wrapper functions at each call.  t is an
 * exception-- it is set by the user just once. levmar calls
 * DFP_check() (below) which either sets t, or points the
 * DFP *dat directly to t, if the wrapper functions are not
 * used (ie for pure C fit functions).
 */

// Don't remember why i called it DFP
typedef struct {
  SV * p_sv; // sv to pass directly to perl fit subs
  SV * d_sv;
  SV * x_sv;
  SV * t_sv;
  pdl * p_pdl; // the pdl pointer
  pdl * d_pdl;
  pdl * x_pdl;
  pdl * t_pdl;
  SV* perl_fit_func; // perl refs to user's fit and jac functions
  SV* perl_jac_func;
  int datatype;
} DFP;

/*  moved to xs
// called from  sub levmar
DFP * _DFP_create () {
  DFP * dat;
  dat = (DFP *)malloc(sizeof( DFP ));
  if ( NULL == dat ) {
    fprintf(stderr, "Can't allocate storage for dat in DFP_create\n");
    exit(1);
  }
  return dat;
}
*/

/*
// called from sub levmar
int _DFP_free (DFP *dat) {
  if ( dat ) {
    free(dat);
  }
}
*/

/* moved into xs
void DFP_set_perl_funcs ( DFP *dat, SV* perl_fit_func, SV* perl_jac_func ) {
  if ( dat == NULL ) {
    fprintf(stderr, "DFP_set_perl_funcs got null struct\n");
    exit(1);
  }
  dat->perl_fit_func = perl_fit_func;
  dat->perl_jac_func = perl_jac_func;
}
*/

void DFP_create_pdls( DFP *dat, int data_type, int m, int n, int nt) {

  pdl *p_pdl, *x_pdl, *d_pdl, *t_pdl;
  SV  *p_sv, *x_sv, *d_sv, *t_sv;
  
  PDL_Long mf_dims[] = {m};
  int num_mf_dims = 1;

  PDL_Long mjac_dims[] = {m,n}; // for jacobian 'd' variable
  int num_mjac_dims = 2;

  PDL_Long n_dims[] = {n};
  int num_n_dims = 1;

  PDL_Long nt_dims[] = {nt};
  int num_nt_dims = 1;


  // create pdls, but with no data;
  p_pdl = pdl_wrap(NULL, data_type, mf_dims, num_mf_dims, delete_levmar_pdls, 0);
  d_pdl = pdl_wrap(NULL, data_type, mjac_dims, num_mjac_dims, delete_levmar_pdls, 0);
  x_pdl = pdl_wrap(NULL, data_type, n_dims, num_n_dims, delete_levmar_pdls, 0);
  t_pdl = pdl_wrap(NULL, data_type, nt_dims, num_nt_dims, delete_levmar_pdls, 0);
  
  p_sv = getref_pdl(p_pdl);
  d_sv = getref_pdl(d_pdl);
  x_sv = getref_pdl(x_pdl);
  t_sv = getref_pdl(t_pdl);

  // otherwise we get a memory leak
  sv_2mortal(p_sv);
  sv_2mortal(d_sv);
  sv_2mortal(x_sv);
  sv_2mortal(t_sv);

  dat->p_pdl = p_pdl;
  dat->d_pdl = d_pdl;
  dat->t_pdl = t_pdl;
  dat->x_pdl = x_pdl;

  dat->p_sv = p_sv;
  dat->d_sv = d_sv;
  dat->t_sv = t_sv;
  dat->x_sv = x_sv;
  
}

/* If dat is null then no DFP struct was created because we
 * are using a C function and it will expect t as the last
 * argument (and dat is passed as the last argument).  If
 * dat is non-null, it will hold a struct with all the pdls
 * and we must put t in its proper place.
 */
void DFP_check(DFP **dat, int data_type, int m, int n, int nt, void *t ) {
  if ( *dat != NULL) {
    DFP_create_pdls( *dat, data_type, m, n, nt);
    (*dat)->t_pdl->data = t;
  }
  else { // Pure C routine that don't need this wrapper, expect t as last arg
    *dat = (DFP *)t; // cast in case compiler complains.
  }
}


// Code from  PDL::API docs. For creating a pdl and giving it
// data storage allocated elswhere.
static pdl* pdl_wrap(void *data, int datatype, PDL_Long dims[],
		     int ndims, DelMagic delete_magic, int delparam)
{
  pdl* npdl = PDL->pdlnew(); /* get the empty container */
  PDL->setdims(npdl,dims,ndims); /* set dims      */
  npdl->datatype = datatype;     /* and data type */
  npdl->data = data;             /* point it to your data */
  /* make sure the core doesn't meddle with your data */
  npdl->state |= PDL_DONTTOUCHDATA | PDL_ALLOCATED;
     if (delete_magic != NULL) {
       PDL->add_deletedata_magic(npdl, delete_magic, delparam);
     }
     else
       PDL->add_deletedata_magic(npdl, default_magic, 0);
     return npdl;
}

// Don't free data, it is allocated inside liblevmar
static void delete_levmar_pdls(pdl* p, int param)
{
  if (p->data) {
    // free(p->data); 
  }
  else {
     }
  p->data = 0;
}

/* Had to cut and paste from pdlcore.c
 * Is this routine somehow otherwise available ?
 * This returns  SV *ref, giving access to the
 * pdl *it as a perl scalar. 
 */
static SV *getref_pdl(pdl *it) {
	SV *newref;
	if(!it->sv) {
		SV *ref;
		HV *stash = gv_stashpv("PDL",TRUE);
		SV *psv = newSViv(PTR2IV(it));
		it->sv = psv;
		newref = newRV_noinc(it->sv);
		(void)sv_bless(newref,stash);
	} else {
		newref = newRV_inc(it->sv);
		SvAMAGIC_on(newref);
	}
	return newref;
}


//-----------------------------------------------------
// Modified Function from pdl gsl interp code
// It follows the perlembed doc example closely

void LEVFUNC(double *p, double *x, int m, int n, DFP *dat) {
  int count;
  dSP;  

  dat->p_pdl->data = p; 
  dat->x_pdl->data = x; 

  ENTER;
  SAVETMPS;

  PUSHMARK(SP);

  // Dont make mortal, they come from outside this routine
  XPUSHs(dat->p_sv);
  XPUSHs(dat->x_sv);
  XPUSHs(dat->t_sv);

  PUTBACK;
  count=call_sv(dat->perl_fit_func,G_SCALAR);
  SPAGAIN; 
  if (count!=1)
    croak("error calling perl function\n");

  PUTBACK;
  FREETMPS;
  LEAVE;

}

void JLEVFUNC(double *p, double *d, int m, int n, DFP *dat) {
  int count;
  dat->p_pdl->data = p; 
  dat->d_pdl->data = d; 

  dSP;  
  ENTER;
  SAVETMPS;
  PUSHMARK(SP);
  XPUSHs(dat->p_sv);
  XPUSHs(dat->d_sv);
  XPUSHs(dat->t_sv);
  PUTBACK;
  count=call_sv(dat->perl_jac_func,G_SCALAR);
  SPAGAIN; 
  if (count!=1)
    croak("error calling perl function\n");
  PUTBACK;
  FREETMPS;
  LEAVE;
}