%module IntegratorLIB
%{
#include "../io/PDBReader.h"
#include "../io/PDBWriter.h"
#include "../frontend/buildTopology.h"
#include "/afs/nd.edu/user25/tcickovs/Research/mdl/internal/extralibs/NUMPY/lib/python2.4/site-packages/numpy/core/include/numpy/ndarrayobject.h"
#include "../base/Real.h"
#include "../topology/CubicCellManager.h"
%}

%include "../frontend/buildTopology.i"
%include "../base/Vector3DBlock.i"
%include "../base/ScalarStructure.i"
%include "../topology/Topology.i"
%include "../forces/ForceGroup.i"
%include "../forces/AngleSystemForce.i"
%include "../forces/BondSystemForce.i"
%include "../forces/DihedralSystemForce.i"
%include "../forces/DihedralIIISystemForce.i"
%include "../forces/ImproperSystemForce.i"
%include "../forces/HarmDihedralSystemForce.i"
%include "../forces/NonbondedSimpleFullSystemForce.i"
%include "../forces/NonbondedFullSystemForce.i"
%include "../forces/NonbondedCutoffSystemForce.i"
%include "../base/PDB.i"
%include "../forces/NonbondedFullEwaldSystemForce.i"
%include "../forces/NonbondedPMEwaldSystemForce.i"
%include "../forces/NonbondedMultiGridSystemForce.i"
%include "../forces/WrapperMetaForce.i"
%include "../base/mathutilities.i"
%include "../io/PDBReader.i"
%include "../io/PDBWriter.i"

%extend std::basic_string<char> {
   PyObject val() {
      return PyObject(self);
   }
}
%extend ProtoMol::PDB {
   ProtoMol::PDBAtom getAtom(int i) {return (self->atoms[i]);}
}
%extend ProtoMol::PDBReader {
   ProtoMol::PDBAtom getAtom(int i) {return self->getAtoms()[i];}
}
%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;
   }

   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;
   }
   
};
%extend ProtoMol::Vector3D{
   Vector3D& copy(const Vector3D& rhs) {
      self->x = rhs.x;
      self->y = rhs.y;
      self->z = rhs.z;
      return (*self);
   }
   void val() {
      printf("[%f %f %f]\n", self->x, self->y, self->z);
   }
   void enter() {
      scanf("%lf %lf %lf", &(self->x), &(self->y), &(self->z));
   }
   Real __getitem__(int oprnd2){
      return (*self)[oprnd2];
   }
   void __setitem__(int oprnd2, Real value) {
     (*self)[oprnd2] = value;
   }
   void setX(Real val) {
      self->x = val;
   }
   void setY(Real val) {
      self->y = val;
   }
   void setZ(Real val) {
      self->z = val;
   }
   Real getX() {
      return self->x;
   }
   Real getY() {
      return self->y;
   }
   Real getZ() {
      return self->z;
   }
   void fromArray(PyObject* rhs) {
      self->x = PyFloat_AsDouble(PySequence_GetItem(rhs, 0));
      self->y = PyFloat_AsDouble(PySequence_GetItem(rhs, 1));
      self->z = PyFloat_AsDouble(PySequence_GetItem(rhs, 2));
   }
    PyObject* toArray() {
        int dims[1] = { 3 };
        float* vec = new float[3];
        vec[0] = self->x;
        vec[1] = self->y;
        vec[2] = self->z;
        PyObject* rhs = PyArray_FromDimsAndData(1,dims,PyArray_FLOAT,(char*)vec);
        return rhs;
    }
   void copyto(PyObject* rhs) {
        PySequence_SetItem(rhs, 0, Py_BuildValue("f", self->x));
        PySequence_SetItem(rhs, 1, Py_BuildValue("f", self->y));
        PySequence_SetItem(rhs, 2, Py_BuildValue("f", self->z));
        if (PySequence_Size(rhs) > 3)
           for (unsigned int i = 3; i < PySequence_Size(rhs); i++)
              PySequence_DelItem(rhs, i);
   }
   Vector3D __add__(Vector3D op1, Vector3D op2) {
      ProtoMol::Vector3D retval;
      retval.x = op1.x + op2.x;
      retval.y = op1.y + op2.y;
      retval.z = op1.z + op2.z;
   }
   Vector3D __sub__(Vector3D op1, Vector3D op2) {
      ProtoMol::Vector3D retval;
      retval.x = op1.x - op2.x;
      retval.y = op1.y - op2.y;
      retval.z = op1.z - op2.z;
   }
 };