/************************************************************************

  Copyright
  Alessandro MIRONE
  mirone@esrf.fr

  Copyright 2002  by European Synchrotron Radiation Facility, Grenoble, 
                  France

                               ----------
 
                           All Rights Reserved
 
                               ----------

Permission to use, copy, modify, and distribute this software and its
documentation for any purpose and without fee is hereby granted,
provided that the above copyright notice appear in all copies and that
both that copyright notice and this permission notice appear in
supporting documentation, and that the names of European Synchrotron
Radiation Facility or ESRF or BLISS not be used in advertising or 
publicity pertaining to distribution of the software without specific, 
written prior permission.

EUROPEAN SYNCHROTRON RADIATION FACILITY DISCLAIMS ALL WARRANTIES WITH
REGARD TO THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF
MERCHANTABILITY AND FITNESS, IN NO EVENT SHALL EUROPEAN SYNCHROTRON
RADIATION FACILITY OR ESRF BE LIABLE FOR ANY SPECIAL, INDIRECT OR 
CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, 
DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER
TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR
PERFORMANCE OF THIS SOFTWARE.

**************************************************************************/
// ************************************
// Alessandro MIRONE 
// Peter      KRUEGER
//
// September 2002
// *************************************


typedef double Hxx_FLOAT;



%module Hilbertxx
%{
      #include "Hilbertxx.hh"
      #include<gsl/gsl_sf_coupling.h>
%}



%{
#include <Numeric/arrayobject.h> 
%}

%include typemaps.i

// THE FOLLOWING TYPEMAPS HAVE BEEN STOLEN FROM
// THE SWIG *.I FILES OF THE PLPLOT PROJECT

%{
   // global variables for consistency check on argument patterns
  static int Alen = 0;
//  static int Xlen = 0, Ylen = 0;
%}



%typemap(in) (int  n, int * Array) (PyArrayObject* tmp) {
  tmp = (PyArrayObject *)PyArray_ContiguousFromObject($input, PyArray_INT, 1, 1);
  if(tmp == NULL) return NULL;
  $1 = Alen = tmp->dimensions[0];
  $2 = (int *)tmp->data;
}
%typemap(freearg) (int n,int * Array) {Py_DECREF(tmp$argnum);}


%typemap(in) ( int * ArrayCHECK) (PyArrayObject* tmp) {
  tmp = (PyArrayObject *)PyArray_ContiguousFromObject($input, PyArray_INT, 1, 1);
  if(tmp == NULL) return NULL;
  if(tmp->dimensions[0] != Alen) {
    PyErr_SetString(PyExc_ValueError, "Vectors must be same length.");
    return NULL;
  }
  $1 = (int *)tmp->data;
}
%typemap(freearg) (int * ArrayCHECK) {Py_DECREF(tmp$argnum);}





%typemap(in) ( int * ArrayNcreation) (PyArrayObject* tmp) {
  tmp = (PyArrayObject *)PyArray_ContiguousFromObject($input, PyArray_INT, 1, 1);
  if(tmp == NULL) return NULL;
  if(tmp->dimensions[0] != (arg1)->Ncreation) {
    PyErr_SetString(PyExc_ValueError, "Vectors for creation operator has not the right lenght");
    return NULL;
  }
  $1 = (int *)tmp->data;
}
%typemap(freearg) (int * ArrayNcreation) {Py_DECREF(tmp$argnum);}



%typemap(in) ( int * ArrayNdestruction) (PyArrayObject* tmp) {
  tmp = (PyArrayObject *)PyArray_ContiguousFromObject($input, PyArray_INT, 1, 1);
  if(tmp == NULL) return NULL;
  if(tmp->dimensions[0] != (arg1)->Ndestruction) {
    PyErr_SetString(PyExc_ValueError, "Vectors for destruction operator has not the right lenght");
    return NULL;
  }
  $1 = (int *)tmp->data;
}
%typemap(freearg) (int * ArrayNdestruction) {Py_DECREF(tmp$argnum);}







%typemap(in) ( int * ArrayNOCHECK) (PyArrayObject* tmp) {
  tmp = (PyArrayObject *)PyArray_ContiguousFromObject($input, PyArray_INT, 1, 1);
  if(tmp == NULL) return NULL;
  if(tmp->dimensions[0] != Alen) {
    PyErr_SetString(PyExc_ValueError, "Vectors must be same length.");
    return NULL;
  }
  $1 = (int *)tmp->data;
}
%typemap(freearg) (int * ArrayNOCHECK) {Py_DECREF(tmp$argnum);}

%init %{	
  import_array();
%}	


%typemap(in) PyObject *pyfunc{
  if (!PyCallable_Check($input)){
    PyErr_SetString(PyExc_TypeError, "Need a callable object!");
    return NULL;
  }
  $1 = $input;
}


class  Hxx_TypeForBits {
public:

  inline void setTo0() ;

  void print_values() ;

  inline void setToOne();
  
  inline void setToValue(int v) ;

  inline int isZero() ;

  int countOne(); 
}  ;


%extend Hxx_TypeForBits {
        char* __str__ () {
	  static char value[sizeof(Hxx_TypeForBits)*8+1];
	  // BE careful this is not thread safe...
	  self->return_string(value);
	  return value;
	}
	    
	Hxx_TypeForBits  __or__ (  Hxx_TypeForBits  &b) {
	  return *self | b ;
	}

	Hxx_TypeForBits  __and__ (  Hxx_TypeForBits  &b) {
	  return *self & b ;
	}
	
	Hxx_TypeForBits  __xor__ (  Hxx_TypeForBits  &b) {
	  return *self ^ b ;
	}

	Hxx_TypeForBits  __rshift__ (  int  b) {
	  return *self >> b ;
	}

	Hxx_TypeForBits  __lshift__ (  int  b) {
	  return *self << b ;
	}

	int  __eq__ ( Hxx_TypeForBits  &b ) {
	  return *self == b ;
	}
	
	int __ne__ ( Hxx_TypeForBits  &b ) {
	  return *self != b ;
	}

	int __lt__ (  Hxx_TypeForBits  &b ) {
	  return *self < b ;
	}

}



/**
 * @short[ Objects of this class represent a second quantisation state ]
 * Objects of this class represent a second quantisation state
 * as a sequence of 1 and 0, like for example
 *     100101011
 * Which are the occupancies. Together with the occupancies
 * the signes to extract a given fermion are given as the integral modulus 2 
 * of the occupancies
 *     011100110
 */
class Hxx_FermionicState {
public:
  
  Hxx_FermionicState(Hxx_TypeForBits occupancies);

  void print_values(int mask=0);

  int  n_particles();

  Hxx_TypeForBits& get_value();

  Hxx_TypeForBits occupancies;
};


class Hxx_Filter {

public:
  void addCondition (Hxx_TypeForBits &mask, char operation, double factor);
  Hxx_Filter (int);
  ~Hxx_Filter ();
  int filtre(Hxx_TypeForBits &state);

  struct opFilter {
    Hxx_TypeForBits mask;
    char operation;
    double factor;
  } opFilter ;

  double min;
  double max;

};




class Hxx_1p_FermionicOperator_  {

public:
  inline Hxx_1p_FermionicOperator_(int pos);
  inline Hxx_TypeForBits & get_1p_state() ;

};

%extend  Hxx_1p_FermionicOperator_ {
  Hxx_1p_FermionicOperator_ * __getitem__ (int pos) {
    
    return self+pos ;
    }
}

class Hxx_1p_FermionicOperators  {

public:
//  inline Hxx_1p_FermionicOperator_ get_operator(int n) ;
  
  static Hxx_1p_FermionicOperator_ * Initialise();
};


class Hxx_N_counter  {
  
public:
  
  inline Hxx_N_counter( int n, int *Array, int reference_n=0) ;

  int N(Hxx_FermionicState* orig);

};



class Hxx_Tree ;



class Hxx_Normal_Operators_Collection {
public:

  Hxx_Normal_Operators_Collection( int Ncreation, int Ndestruction,
				   int chunck_size = 10000);

  void add_product(int * ArrayNcreation /* creators */ , int *ArrayNdestruction /* destructors */, Hxx_FLOAT coeff );
  void reorder();
  
  inline int get_n_items() ;
  
  Hxx_FLOAT operate( int i, Hxx_FermionicState* orig, Hxx_FermionicState* result, int & INOUT  /* next */);
  %name(myprint) void print();

};

// -------------------------------------------------------------------------------------------

%{
  /* overload filter function as a callback. This in order to write it in Python and use it in C++ */

  static int PyFilter (Hxx_FermionicState * s, void *ref)
    {
      PyObject *func, *arglist;
      PyObject *result;
      PyObject *arg1;

      int ires = 0;

      func = (PyObject *) ref;

      arg1 = SWIG_NewPointerObj((Hxx_FermionicState *) s, SWIGTYPE_p_Hxx_FermionicState, 0);
      arglist=Py_BuildValue("(O)",arg1);
      
      Py_DECREF(arg1);
      
      result = PyEval_CallObject(func, arglist);

      Py_DECREF(arglist);
      // Py_DECREF(func);
      
      if (result) {
	ires = PyInt_AsLong(result);
      }
      Py_XDECREF(result);

      return ires;
	}
  %}


typedef int (*Callback) (Hxx_FermionicState *, void *);

class Hxx_basis {
public:
  Hxx_basis(int initialdimension, int chunck_size=10000);
  void add_state(Hxx_FermionicState * state);

  int operate_andputin_Tree( Hxx_Normal_Operators_Collection * collection, int Ncolls=1, Hxx_Filter* filtre=0 );

  void reconstruct_basis();
  void set_filter(Callback filtre, void *clientdata);

  void reconstruct_basis(CounterForHxxBasisFilter *);

  void reconstruct_basis(Hxx_Filter* Hxfiltre);	


  inline int get_dimension() { return this->nstates; }

  int find_state(Hxx_FermionicState * state  );

  inline Hxx_FermionicState *get_basis() { return this->basis;}

private:
  int nstates;
  int Ncreated;
  Hxx_FermionicState *basis ;
  Callback  filtre;
  void *clientdata;
  Hxx_Tree  * orderer ;
  Hxx_Chunker * states_chunker;
  double average_depth;
  double average_depth2;
};

%extend Hxx_basis {
  void set_pyfilter(PyObject *pyfunc){
    self->set_filter(PyFilter, (void *) pyfunc);
    Py_INCREF(pyfunc);
  }
}

// this method to acced states contained in array of fermionic states
// like the one returned by Hxx_basis::get_basis
%extend Hxx_FermionicState {
	Hxx_FermionicState * __getitem__(int j) {
		return self+j;
	}	
};




// -----------------------------------------------------------------------------------------

class Hxx_transition {
public:
  
  inline Hxx_transition(int from , int to, Hxx_FLOAT coeff) ;
  
  inline void initialise(int from , int to, Hxx_FLOAT coeff) ;
  
  inline void add_coeff(Hxx_FLOAT coeff) ;
  
  // static  int less(Hxx_transition *a, Hxx_transition *b) ;

  inline int get_from() ;
  inline int get_to() ;
  inline Hxx_FLOAT  get_coeff() ;
  

};



%typemap(in) (int n_counters, Hxx_N_counter  ** diag_pro)  {
  $1 = PyList_Size( $input );
  $2 = new Hxx_N_counter  * [ $1 ] ;
  for(int i=0; i< $1 ; i++ ) {
    SWIG_ConvertPtr(PyList_GetItem( $input ,  i ) ,(void**) &($2[i]),
		    SWIGTYPE_p_Hxx_N_counter ,0);
  }  
}

%typemap(freearg) (int n_counters, Hxx_N_counter  ** diag_pro) {
  free($2);
}



%typemap(in, numinputs=0) ( int  *& from,  int  *& to ,Hxx_FLOAT *& coeff, int & ntrans ) (int *ia, int *ib, Hxx_FLOAT * coeF, int N_tmp) {
  $1 = &ia;
  $2 = &ib;
  $3 = &coeF ;
  $4 = &N_tmp;
} 
%typemap(argout) ( int  *& from,  int  *& to , Hxx_FLOAT *& coeff,int & ntrans  )  {
  int nd=1;
  int dims[1];
  PyObject *arr_a, *arr_b, *arr_c;
  dims[0]=N_tmp$argnum;

  arr_a = PyArray_FromDims(nd,dims,PyArray_INT) ;
  memcpy( ( (PyArrayObject *) arr_a )->data  ,  *$1 , dims[0]*sizeof(int) );
  arr_b = PyArray_FromDims(nd,dims,PyArray_INT) ;
  memcpy( ( (PyArrayObject *) arr_b )->data  ,  *$2 , dims[0]*sizeof(int) );
  arr_c = PyArray_FromDims(nd,dims,PyArray_DOUBLE) ;
  memcpy( ( (PyArrayObject *) arr_c )->data  ,  *$3 , dims[0]*sizeof(double) );

  delete [] *$1;
  delete [] *$2;
  delete [] *$3;

  resultobj = t_output_helper(resultobj,arr_a );
  resultobj = t_output_helper(resultobj,arr_b );
  resultobj = t_output_helper(resultobj,arr_c );

} 



class Hxx_TransitionMatrix {
public:

  Hxx_TransitionMatrix();

  void add_contribution(Hxx_basis * a_basis, Hxx_basis *b_basis ,
			Hxx_Normal_Operators_Collection  * Norm_Op, Hxx_FLOAT coeff=1.0, int complete_basis=1);

  void add_contribution_debug(Hxx_basis * a_basis, Hxx_basis *b_basis ,
			Hxx_Normal_Operators_Collection  * Norm_Op, Hxx_FLOAT coeff=1.0, int complete_basis=1);

  void add_contribution(Hxx_basis * a_basis, Hxx_basis *b_basis ,
			int n_counters, Hxx_N_counter  ** diag_pro, Hxx_FLOAT  coeff);  // <<<<<< 

  void Write_and_Clean(char *nome_file, int binary=1);
  void Read_add(char *nome_file, Hxx_FLOAT, int binary=1);

  void create_3Array_and_Clean(int *& from, int *& to, Hxx_FLOAT *& coeff, int & ntrans  );  // <<<<< 
  
  void free_space();
  
};





%name(Wigner3J) double gsl_sf_coupling_3j( int l1, int l2 , int l3 ,   int m1 , int m2 , int m3    );
