%module DMDLeapfrogIntegrator
%{
#include "DMDLeapfrogIntegrator.h"
#include "../base/Vector3DBlock.h"
%}

%rename (ProtoMol_IntegratorDefinition_print) ProtoMol::IntegratorDefinition::print;

%include "std_vector.i"
%template() std::vector<ProtoMol::Vector3D >;
%include "../base/Real.h"
%include "../base/Proxy.h"
%include "../base/simpleTypes.h"
%include "../base/Vector3D.h"
%include "../base/Vector3DBlock.h"
%include "../base/AbstractEnumType.h"
%include "../base/ConstraintValueType.h"
%include "../base/ValueType.h"
%include "../base/Value.h"
%include "../base/Parameter.h"
%include "../base/MakeableDefinition.h"
%include "../base/Makeable.h"
%include "../base/XSC.h"

%include "IntegratorDefinition.h"
%include "Integrator.h"
%include "StandardIntegrator.h"
%include "STSIntegrator.h"
%include "DMDLeapfrogIntegrator.h"


%extend ProtoMol::Vector3DBlock {
   void copyto(PyObject* rhs) {
      for (unsigned int i = 0; i < self->size(); i++) {
        PySequence_SetItem(rhs, 3*i+0, Py_BuildValue("f", self->myData[3*i]));
        PySequence_SetItem(rhs, 3*i+1, Py_BuildValue("f", self->myData[3*i+1]));
        PySequence_SetItem(rhs, 3*i+2, Py_BuildValue("f", self->myData[3*i+2]));
      }
      if (PySequence_Size(rhs) > 3*self->size()+3)
        for (unsigned int j = 3*self->size()+3; j < PySequence_Size(rhs); j++)
          PySequence_DelItem(rhs, j);
   }

   void copytoNative(PyObject* rhs) {
      for (unsigned int i = 0; i < self->size(); i++) {
        PySequence_SetItem(PyList_GetItem(rhs, i), 0, Py_BuildValue("f", self->myData[3*i]));
        PySequence_SetItem(PyList_GetItem(rhs, i), 1, Py_BuildValue("f", self->myData[3*i+1]));
        PySequence_SetItem(PyList_GetItem(rhs, i), 2, Py_BuildValue("f", self->myData[3*i+2]));
        if (PySequence_Size(PyList_GetItem(rhs, i)) > 3)
           for (unsigned int j = 3; j < PySequence_Size(PyList_GetItem(rhs, i)); j++)
              PySequence_DelItem(PyList_GetItem(rhs, i), j);
      }  
   }
   void pushback(ProtoMol::Vector3D* vd) {
       self->push_back(*vd);
   }
   void pushback(Real x, Real y, Real z) {
       self->push_back(ProtoMol::Vector3D(x,y,z));
   }
   int length() {
       return self->size();
   }
   void rsize(int a) {
       self->resize(a);
   }
   ProtoMol::Vector3D __getitem__(int oprnd2){
	return ProtoMol::Vector3D(self->myData[oprnd2*3], self->myData[oprnd2*3+1], self->myData[oprnd2*3+2]);	
	//return Py_BuildValue("f", self->myData[oprnd2]);
   }

   void __setitem__(int oprnd2, Real val) {
	self->myData[oprnd2] = val;
   }

   void setData(PyObject* rhs) {
      self->loadArr();
      self->mySize = ((PyArrayObject*)rhs)->dimensions[0]/3;
      self->myArraySize = ((PyArrayObject*)rhs)->dimensions[0];
	if (!(PyArrayObject*)(self->myData)) {
	delete self->myData;
	}
      self->myData = (double*)(((PyArrayObject*)rhs)->data);
   }

   PyObject* getData() {
	int mySize = self->size()*3;
      int dims[1] = {mySize};
	self->loadArr();
	PyObject* rhs = PyArray_SimpleNewFromData(1,dims,PyArray_DOUBLE,(char*)(self->myData));
        return rhs;
   }

   ProtoMol::Vector3D get(int oprnd2){
	return ProtoMol::Vector3D(self->myData[oprnd2*3], self->myData[oprnd2*3+1], self->myData[oprnd2*3+2]);
//      return (*self)[oprnd2];
   }

   void __setitem__(int oprnd2, ProtoMol::Vector3D vec) {
	self->myData[oprnd2*3] = vec.x;
	self->myData[oprnd2*3+1] = vec.y;
	self->myData[oprnd2*3+2] = vec.z;
      //(*self)[oprnd2] = vec;
   }

   Vector3DBlock operator+(Vector3DBlock v1) {
      ProtoMol::Vector3DBlock* retval = new ProtoMol::Vector3DBlock();
      retval->resize(v1.size());
      for (unsigned int i = 0; i < v1.size()*3; i++) {
	retval->myData[i] = v1.myData[i] + self->myData[i];
      }
      return *retval;
   }
   
   Vector3DBlock operator-(Vector3DBlock v1) {
      ProtoMol::Vector3DBlock* retval = new ProtoMol::Vector3DBlock();
      retval->resize(v1.size());
      for (unsigned int i = 0; i < v1.size()*3; i++) {
	retval->myData[i] = self->myData[i] - v1.myData[i];
      }
      return *retval;
   }

   Vector3DBlock operator*(float scale) {
      ProtoMol::Vector3DBlock* retval = new ProtoMol::Vector3DBlock();
      retval->resize(self->size());
      for (unsigned int i = 0; i < self->size(); i++) {
	retval->myData[i] = self->myData[i] * scale;
      }
      return *retval;
   }

   Vector3DBlock operator*(PyObject* rhs) {
      ProtoMol::Vector3DBlock* retval = new ProtoMol::Vector3DBlock();
      retval->resize(self->size());
      if (PyFloat_Check(rhs)){
         for (unsigned int i = 0; i < self->size(); i++) {
		retval->myData[i] = self->myData[i] * PyFloat_AsDouble(rhs);
         }
      }
      else{
      int diagonal = !(PyObject_Not(PyObject_GetAttrString(rhs, "diagonal")));
      if (!diagonal) {
      for (unsigned int i = 0; i < 3; i++) {
        for (unsigned int j = 0; j < self->size(); j++) {
	   for (unsigned int k = 0; k < self->size(); k++) {
	      if (i == 0) (*retval).myData[j*3] += (*self).myData[k*3] * PyFloat_AsDouble(PySequence_GetItem(PySequence_GetItem(rhs, k), j));
	      else if (i == 1) (*retval).myData[j*3+1] += (*self).myData[k*3+1] * PyFloat_AsDouble(PySequence_GetItem(PySequence_GetItem(rhs, k), j));
	      else (*retval).myData[j*3+2] += (*self).myData[k*3+2] * PyFloat_AsDouble(PySequence_GetItem(PySequence_GetItem(rhs, k), j));
	   }
	}
      }
      }
      else  // Optimized version for diagonal matrices
      {
         for (unsigned int i = 0; i < 3; i++) {
	    for (unsigned int j = 0; j < self->size(); j++) {
	      if (i == 0) (*retval).myData[j*3] += (*self).myData[j*3] * PyFloat_AsDouble(PySequence_GetItem(PySequence_GetItem(rhs, j), j));
	      else if (i == 1) (*retval).myData[j*3+1] += (*self).myData[j*3+1] * PyFloat_AsDouble(PySequence_GetItem(PySequence_GetItem(rhs, j), j));
	      else (*retval).myData[j*3+2] += (*self).myData[j*3+2] * PyFloat_AsDouble(PySequence_GetItem(PySequence_GetItem(rhs, j), j));
	    }
	 }
      }
      }
      return *retval;
   }
   
};