/** Makes sure a Numpy Array has the given type and dimension.  Raises
a Python exception if it does not.
@param type_num See http://docs.scipy.org/doc/numpy/reference/c-api.dtype.html eg: NPY_DOUBLE */
PyArrayObject *check_dimensions(
PyObject *ovec,
std::string const &vname,
int type_num,
int ndim,
int const *dims,
int N_ndim)
{
	// Check user and template #-dimensions match
	if (ndim != N_ndim) {
 		char buf[200 + vname.size()];
		sprintf(buf, "check_dimensions(%s): Number of template dimensions (%d) must match number of array dimensions (%d)", vname.c_str(), N_ndim, ndim);
		fprintf(stderr, "%s\n", buf);
		PyErr_SetString(PyExc_ValueError, buf);
		throw std::exception();
	}

	// Check that it's not null
	if (!ovec) {
		std::string serr = "check_dimensions: Array object " + vname + " is null";
		fprintf(stderr, "%s\n", serr.c_str());
		PyErr_SetString(PyExc_ValueError, serr.c_str());
		throw std::exception();
	}

	// Check that it's type PyArrayObject
	if (!PyArray_Check(ovec)) {
		std::string serr = "check_dimensions: Object " + vname + " is not a Numpy array";
		fprintf(stderr, "%s\n", serr.c_str());
		PyErr_SetString(PyExc_ValueError,
			serr.c_str());
		throw std::exception();
	}
	PyArrayObject *vec = (PyArrayObject *)ovec;

	// Check the data type and number of dimensions
	if (vec->descr->type_num != type_num || vec->nd != ndim)  {
		char buf[200 + vname.size()];
		sprintf(buf, "check_dimensions: %s must be of type_num %d and %d dimensions (its is of type_num=%d and %d dimensions).", vname.c_str(), type_num, ndim, vec->descr->type_num, vec->nd);
		fprintf(stderr, "%s\n", buf);
		PyErr_SetString(PyExc_ValueError, buf);
		throw std::exception();
	}

	// Check the dimensions themselves
	for (int i=0; i<ndim; ++i) {
		if (dims[i] < 0) continue;		// Don't check this dimension
		if (dims[i] != vec->dimensions[i]) {
			char buf[200 + vname.size()];
			sprintf(buf,
				"%s: Array dimension #%d is %d, should be %d",
				vname.c_str(), i, vec->dimensions[i], dims[i]);
			fprintf(stderr, "%s\n", buf);
			PyErr_SetString(PyExc_ValueError, buf);
			throw std::exception();
		}
	}
	return vec;
}

static int PyLWPR_S_init_alpha(PyLWPR *self,PyObject *value, void *closure) { 
   if (!PyArray_Check(value)) {
      PyErr_SetString(PyExc_TypeError, "Attribute 'init_alpha' must be a numpy array.");
      return -1;
   }
   return set_matrix_from_array(self->model.nIn, self->model.nInStore, self->model.nIn, 
                                 self->model.init_alpha, (PyArrayObject *)value);
}

    static void convert(PyObject *src, T &dst) {
      if (PyArray_Check(src)) {
	Converter<PyArrayObject*, T>::convert(src, dst);
      }
      else {
	throw std::string("Cannot convert object!");
      }
    }

    static PyArrayObject *cast(PyObject *obj) {
      if (PyArray_Check(obj)) {
	return (PyArrayObject*)obj;
      }
      else {
	throw std::string("Object must be a Python array!");
      }
    }

td_tag_t td_py_get_eltype(void *v)
{
    if (PyArray_Check( (PyObject *) v)){
        PyArrayObject *arr = (PyArrayObject *)v;
        return numpy_type_to_td(PyArray_TYPE(v));
    }
    return TD_UNKNOWN;
}

NPY_NO_EXPORT npy_bool
_IsWriteable(PyArrayObject *ap)
{
    PyObject *base=PyArray_BASE(ap);
#if defined(NPY_PY3K)
    Py_buffer view;
#else
    void *dummy;
    Py_ssize_t n;
#endif

    /* If we own our own data, then no-problem */
    if ((base == NULL) || (PyArray_FLAGS(ap) & NPY_ARRAY_OWNDATA)) {
        return NPY_TRUE;
    }
    /*
     * Get to the final base object
     * If it is a writeable array, then return TRUE
     * If we can find an array object
     * or a writeable buffer object as the final base object
     * or a string object (for pickling support memory savings).
     * - this last could be removed if a proper pickleable
     * buffer was added to Python.
     *
     * MW: I think it would better to disallow switching from READONLY
     *     to WRITEABLE like this...
     */

    while(PyArray_Check(base)) {
        if (PyArray_CHKFLAGS((PyArrayObject *)base, NPY_ARRAY_OWNDATA)) {
            return (npy_bool) (PyArray_ISWRITEABLE((PyArrayObject *)base));
        }
        base = PyArray_BASE((PyArrayObject *)base);
    }

    /*
     * here so pickle support works seamlessly
     * and unpickled array can be set and reset writeable
     * -- could be abused --
     */
    if (PyString_Check(base)) {
        return NPY_TRUE;
    }
#if defined(NPY_PY3K)
    if (PyObject_GetBuffer(base, &view, PyBUF_WRITABLE|PyBUF_SIMPLE) < 0) {
        PyErr_Clear();
        return NPY_FALSE;
    }
    PyBuffer_Release(&view);
#else
    if (PyObject_AsWriteBuffer(base, &dummy, &n) < 0) {
        PyErr_Clear();
        return NPY_FALSE;
    }
#endif
    return NPY_TRUE;
}

/* optimize float array or complex array to a scalar power */
static PyObject *
fast_scalar_power(PyArrayObject *a1, PyObject *o2, int inplace)
{
    double exp;

    if (PyArray_Check(a1) && array_power_is_scalar(o2, &exp)) {
        PyObject *fastop = NULL;
        if (PyArray_ISFLOAT(a1) || PyArray_ISCOMPLEX(a1)) {
            if (exp == 1.0) {
                /* we have to do this one special, as the
                   "copy" method of array objects isn't set
                   up early enough to be added
                   by PyArray_SetNumericOps.
                */
                if (inplace) {
                    Py_INCREF(a1);
                    return (PyObject *)a1;
                } else {
                    return PyArray_Copy(a1);
                }
            }
            else if (exp == -1.0) {
                fastop = n_ops.reciprocal;
            }
            else if (exp ==  0.0) {
                fastop = n_ops.ones_like;
            }
            else if (exp ==  0.5) {
                fastop = n_ops.sqrt;
            }
            else if (exp ==  2.0) {
                fastop = n_ops.square;
            }
            else {
                return NULL;
            }

            if (inplace) {
                return PyArray_GenericInplaceUnaryFunction(a1, fastop);
            } else {
                return PyArray_GenericUnaryFunction(a1, fastop);
            }
        }
        else if (exp==2.0) {
            fastop = n_ops.multiply;
            if (inplace) {
                return PyArray_GenericInplaceBinaryFunction
                    (a1, (PyObject *)a1, fastop);
            }
            else {
                return PyArray_GenericBinaryFunction
                    (a1, (PyObject *)a1, fastop);
            }
        }
    }
    return NULL;
}

 /**
 * Function to call to evaluate model
 * @param args: input numpy array  [x[],y[]]
 * @return: numpy array object 
 */
 static PyObject * evaluateTwoDimXY( SCCrystalModel* model, 
                              PyArrayObject *x, PyArrayObject *y)
 {
    PyArrayObject *result;
    int x_len, y_len, dims[1];
    //check validity of input vectors
    if (x->nd != 1 || x->descr->type_num != PyArray_DOUBLE
        || y->nd != 1 || y->descr->type_num != PyArray_DOUBLE
        || y->dimensions[0] != x->dimensions[0]){
        const char * message= "evaluateTwoDimXY  expect 2 numpy arrays";
        PyErr_SetString(PyExc_ValueError , message); 
        return NULL;
    }
   
	if (PyArray_Check(x) && PyArray_Check(y)) {
		
	    x_len = dims[0]= x->dimensions[0];
        y_len = dims[0]= y->dimensions[0];
	    
	    // Make a new double matrix of same dims
        result=(PyArrayObject *) PyArray_FromDims(1,dims,NPY_DOUBLE);
        if (result == NULL){
	    const char * message= "Could not create result ";
        PyErr_SetString(PyExc_RuntimeError , message);
	    return NULL;
	    }
       
        /* Do the calculation. */
#pragma omp parallel for
        for (int i=0; i< x_len; i++) {
            double x_value = *(double *)(x->data + i*x->strides[0]);
  		    double y_value = *(double *)(y->data + i*y->strides[0]);
  			double *result_value = (double *)(result->data +
  			      i*result->strides[0]);
  			*result_value = (*model)(x_value, y_value);
        }           
        return PyArray_Return(result); 
        
        }else{
		    PyErr_SetString(CSCCrystalModelError, 
                   "CSCCrystalModel.evaluateTwoDimXY couldn't run.");
	        return NULL;
		}      	
}

/*
 * new_points
 *
 * Stores the items of a Numerical Python array or a Python sequence in a
 * C vector of double values.  Returns the number of elements or -1.
 * The C vector must be freed with g_free().
 */
static gint
new_points(PyObject *sequence, gdouble **out_points, int *out_array_type)
{
    int n = 0;
    gdouble *points = NULL;

    *out_array_type = PyArray_NOTYPE;
#ifndef WITHOUT_NUMPY
    if (PyArray_Check(sequence)) {
        PyArrayObject *array;

        array = (PyArrayObject *)
                PyArray_ContiguousFromObject(sequence, PyArray_DOUBLE, 1, 1);
        if (!array)
            return -1;
        n = array->dimensions[0];
        points = g_new(gdouble, n);
        memcpy(points, array->data, n * sizeof(double));
        Py_DECREF(array);
        *out_array_type = ((PyArrayObject *) sequence)->descr->type_num;
    } else
#endif
        if (PySequence_Check(sequence)) {
            n = PySequence_Length(sequence);
            if (n > 0) {
                int i;

                points = g_new(gdouble, n);
                for (i = 0; i < n; ++i) {
                    PyObject *item, *value;

                    item = PySequence_GetItem(sequence, i);
                    if (PyFloat_Check(item)) {
                        points[i] = PyFloat_AS_DOUBLE(item);
                    } else if (PyNumber_Check(item)
                               && (value = PyNumber_Float(item))) {
                        points[i] = PyFloat_AS_DOUBLE(value);
                        Py_DECREF(value);
                    } else {
                        PyErr_SetString(PyExc_TypeError,
                                        "sequence items must be numbers");
                        Py_DECREF(item);
                        g_free(points);
                        points = NULL;
                        return -1;
                    }
                    Py_DECREF(item);
                }
            }
        } else if (sequence != Py_None) {
            PyErr_SetString(PyExc_TypeError, "argument must be sequence or None");
            return -1;
        }
    *out_points = points;
    return n;
}

static int try_pyarr_from_string(PyObject *obj,const string str) {
  PyArrayObject *arr = NULL;
  if (PyArray_Check(obj) && (!((arr = (PyArrayObject *)obj) == NULL)))
    { STRINGCOPYN(arr->data,str,PyArray_NBYTES(arr)); }
  return 1;
capi_fail:
  PRINTPYOBJERR(obj);
  PyErr_SetString(_lbfgsb_error,"try_pyarr_from_string failed");
  return 0;
}

/*NUMPY_API
  Compute the size of an array (in number of items)
*/
NPY_NO_EXPORT intp
PyArray_Size(PyObject *op)
{
    if (PyArray_Check(op)) {
        return PyArray_SIZE((PyArrayObject *)op);
    }
    else {
        return 0;
    }
}

/*NUMPY_API*/
NPY_NO_EXPORT int
PyArray_CopyObject(PyArrayObject *dest, PyObject *src_object)
{
    PyArrayObject *src;
    PyObject *r;
    int ret;

    /*
     * Special code to mimic Numeric behavior for
     * character arrays.
     */
    if (dest->descr->type == PyArray_CHARLTR && dest->nd > 0 \
        && PyString_Check(src_object)) {
        intp n_new, n_old;
        char *new_string;
        PyObject *tmp;

        n_new = dest->dimensions[dest->nd-1];
        n_old = PyString_Size(src_object);
        if (n_new > n_old) {
            new_string = (char *)malloc(n_new);
            memmove(new_string, PyString_AS_STRING(src_object), n_old);
            memset(new_string + n_old, ' ', n_new - n_old);
            tmp = PyString_FromStringAndSize(new_string, n_new);
            free(new_string);
            src_object = tmp;
        }
    }

    if (PyArray_Check(src_object)) {
        src = (PyArrayObject *)src_object;
        Py_INCREF(src);
    }
    else if (!PyArray_IsScalar(src_object, Generic) &&
             PyArray_HasArrayInterface(src_object, r)) {
        src = (PyArrayObject *)r;
    }
    else {
        PyArray_Descr* dtype;
        dtype = dest->descr;
        Py_INCREF(dtype);
        src = (PyArrayObject *)PyArray_FromAny(src_object, dtype, 0,
                                               dest->nd,
                                               FORTRAN_IF(dest),
                                               NULL);
    }
    if (src == NULL) {
        return -1;
    }

    ret = PyArray_MoveInto(dest, src);
    Py_DECREF(src);
    return ret;
}

static td_tag_t py_type_to_td(PyObject *pVal)
{
    if (PyInt_Check(pVal)){
#if INT32_MAX==LONG_MAX
        return TD_INT32;
#else
        return longlong_to_td(PyInt_AS_LONG(pVal));
#endif
    }
    else if (PyLong_Check(pVal)){
        PyObject *pErr;
        long long v;
        v = PyLong_AsLongLong(pVal);
        if ( v == -1 && PyErr_Occurred()) {
            // could not eval as long long try unsigned long long later
            PyErr_Clear();
        }
        else {
        	return longlong_to_td(PyInt_AS_LONG(v));
        }

        unsigned long long vv;
        vv = PyLong_AsUnsignedLongLong(pVal);
        if ( v == (unsigned long long) -1 && PyErr_Occurred() ) {
            PyErr_Clear();
            return TD_OBJECT;
        }
        else {
            return TD_UINT64;
        }
    }
    else if (PyFloat_Check(pVal)) {
        return TD_DOUBLE;
    }
    else if (PyUnicode_Check(pVal)){
        PyObject *pStr;
        pStr = PyUnicode_AsUTF8String(pVal);
        if ( pStr && PyErr_Occurred() ){
            PyErr_Clear();
            return TD_OBJECT;
        }
        else {
            return TD_UTF8;
        }
    }
    else if (PyString_Check(pVal)) {
        return TD_UTF8;
    }
    else if (PyArray_Check(pVal)) {
        return TD_ARRAY;
    }
    return TD_OBJECT;
}

PyArrayObject *call_funGetData(PyObject *func, double size) 
{
	PyObject *args;
	PyObject *kwargs;
	PyObject *result;
	PyArrayObject *retval;

	/* Make sure we own the GIL */
	PyGILState_STATE state = PyGILState_Ensure();

	/* Verify that func is a proper callable */
	if (!PyCallable_Check(func)) {
		fprintf(stderr, "call_func: expected a callable\n");
		goto fail;
	}
	/* Build arguments */
	args = Py_BuildValue("(d)", size);
	kwargs = NULL;

	/* Call the function */
	result = PyObject_Call(func, args, kwargs);
	Py_DECREF(args);
	Py_XDECREF(kwargs);

	/* Check for Python exceptions (if any) */
	if (PyErr_Occurred()) {
		PyErr_Print();
		goto fail;
	}

  /* Verify the result is a float object */
  if (!PyArray_Check( result )) {
    fprintf( stderr, "call_func: callable didn't return a array\n" );
    goto fail;
  }

  retval = (PyArrayObject*)PyArray_FROM_OTF( result, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY );

	/* Create the return value */
	//retval = (PyArrayObject *)(result);
  //col = (double*)PyArray_GETPTR2( retval, TICA_FIRST_ROW, 0);
	Py_DECREF(result);

	/* Restore previous GIL state and return */
	PyGILState_Release(state);
	return retval;

fail:
	Py_XDECREF(result);
	PyGILState_Release(state);
	return NULL;
	abort();   // Change to something more appropriate
}

// assume numDerivVars x numDerivVars as helper in Hessian conversion
// and lower triangular storage in Hessians
// returns false if conversion failed
bool NRELPythonApplicInterface::python_convert(PyObject *pym, 
					   RealSymMatrix &rm)
{
  // for now, the numpy case isn't called (since handled in calling
  // Hessian array convert)
#ifdef DAKOTA_PYTHON_NUMPY
  if (userNumpyFlag) {
    if (!PyArray_Check(pym) || PyArray_NDIM(pym) != 2 || 
	PyArray_DIM(pym,0) != numDerivVars  ||
	PyArray_DIM(pym,1) != numDerivVars) {
      Cerr << "Python numpy array not 2D of size " << numDerivVars << "x" 
	   << numDerivVars << "." << std::endl;
      return(false);
    }
    PyArrayObject *pao = (PyArrayObject *) pym;
    for (int i=0; i<numDerivVars; ++i)
      for (int j=0; j<=i; ++j)
	rm(i,j) = *(double *)(pao->data + i*(pao->strides[0]) + 
			       j*(pao->strides[1]));
  }
  else
#endif
  {
    if (!PyList_Check(pym) || PyList_Size(pym) != numDerivVars) {
      Cerr << "Python matrix must have " << numDerivVars << "rows." <<std::endl;
      return(false);
    }
    PyObject *pyv, *val;
    for (int i=0; i<numDerivVars; ++i) {
      pyv = PyList_GetItem(pym, i);
      if (!PyList_Check(pyv) || PyList_Size(pyv) != numDerivVars) {
	Cerr << "Python vector must have length " << numDerivVars << "." 
	     << std::endl;
	return(false);
      }
      for (int j=0; j<=i; ++j) {
	val = PyList_GetItem(pyv, j);
	if (PyFloat_Check(val))
	  rm(i,j) = PyFloat_AsDouble(val);
	else if (PyInt_Check(val))
	  rm(i,j) = (double) PyInt_AsLong(val);
	else {
	  Cerr << "Unsupported Python data type converting vector." 
	       << std::endl;
	  Py_DECREF(val);
	  return(false);
	}
      }
    }

  }
  return(true);
}

static PyObject *zero_if_zero_inplace(PyObject *self, PyObject *args)
{
  PyObject* a = NULL;
  PyObject* b = NULL;
  npy_intp sz = 0, i;
  npy_complex64* tmp = NULL;
  npy_complex64* a_data = NULL;
  npy_float32* b_data = NULL;
  if (!PyArg_ParseTuple(args, "OO", &a, &b))
    return NULL;
  if (!(PyArray_Check(a) && PyArray_Check(b)))
    {
      PyErr_SetString(PyExc_TypeError,"arguments must be array objects");
      return NULL;
    }
  sz = PyArray_SIZE(a);

  if (sz != PyArray_SIZE(b))
    {
      PyErr_SetString(PyExc_TypeError,"argument sizes must be equal");
      return NULL;
    }
  if (! ((PyArray_TYPE(a) == PyArray_COMPLEX64) && (PyArray_TYPE(b) == PyArray_FLOAT32)))
    {
      PyErr_SetString(PyExc_TypeError,"argument types must be complex64 and float32");
      return NULL;
    }
  a_data = (npy_complex64*)PyArray_DATA(a);
  b_data = (npy_float32*)PyArray_DATA(b);
  for (i=0; i<sz; ++i)
    {
      if (b_data[i]==0)
	{
	  tmp = a_data + i;
	  tmp->real = tmp->imag = 0.0;
	}
    }
  return Py_BuildValue("");
}

fvec_t *
PyAubio_ArrayToCFvec (PyObject *input) {
  PyObject *array;
  fvec_t *vec;
  if (input == NULL) {
    PyErr_SetString (PyExc_ValueError, "input array is not a python object");
    goto fail;
  }
  // parsing input object into a Py_fvec
  if (PyArray_Check(input)) {

    // we got an array, convert it to an fvec 
    if (PyArray_NDIM (input) == 0) {
      PyErr_SetString (PyExc_ValueError, "input array is a scalar");
      goto fail;
    } else if (PyArray_NDIM (input) > 1) {
      PyErr_SetString (PyExc_ValueError,
          "input array has more than one dimensions");
      goto fail;
    }

    if (!PyArray_ISFLOAT (input)) {
      PyErr_SetString (PyExc_ValueError, "input array should be float");
      goto fail;
    } else if (PyArray_TYPE (input) != AUBIO_NPY_SMPL) {
      PyErr_SetString (PyExc_ValueError, "input array should be float32");
      goto fail;
    } else {
      // input data type is float32, nothing else to do
      array = input;
    }

    // vec = new_fvec (vec->length);
    // no need to really allocate fvec, just its struct member 
    vec = (fvec_t *)malloc(sizeof(fvec_t));
    vec->length = PyArray_SIZE (array);
    vec->data = (smpl_t *) PyArray_GETPTR1 (array, 0);

  } else if (PyObject_TypeCheck (input, &PyList_Type)) {
    PyErr_SetString (PyExc_ValueError, "does not convert from list yet");
    return NULL;
  } else {
    PyErr_SetString (PyExc_ValueError, "can only accept vector of float as input");
    return NULL;
  }

  return vec;

fail:
  return NULL;
}

 void transConformer(Conformer &conf, python::object trans) {
   PyObject *transObj = trans.ptr();
   if (!PyArray_Check(transObj)) {
     throw_value_error("Expecting a numeric array for transformation");
   }
   PyArrayObject *transMat = reinterpret_cast<PyArrayObject *>(transObj);
   unsigned int nrows = transMat->dimensions[0];
   unsigned int dSize = nrows*nrows;
   double *inData = reinterpret_cast<double *>(transMat->data);
   RDGeom::Transform3D transform;
   double *tData  = transform.getData();
   memcpy(static_cast<void *>(tData), static_cast<void *>(inData), dSize*sizeof(double));
   MolTransforms::transformConformer(conf, transform);
 }

/* Converts a Python array-like object into a behaved input/output array. */
static int
NI_ObjectToInputOutputArray(PyObject *object, PyArrayObject **array)
{
    /*
     * This is also done in NI_ObjectToOutputArray, double checking here
     * to provide a more specific error message.
     */
    if (PyArray_Check(object) &&
            !PyArray_ISWRITEABLE((PyArrayObject *)object)) {
        PyErr_SetString(PyExc_ValueError, "input/output array is read-only.");
        return 0;
    }
    return NI_ObjectToOutputArray(object, array);
}

static PyObject *
arrayflags_new(PyTypeObject *NPY_UNUSED(self), PyObject *args, PyObject *NPY_UNUSED(kwds))
{
    PyObject *arg=NULL;
    if (!PyArg_UnpackTuple(args, "flagsobj", 0, 1, &arg)) {
        return NULL;
    }
    if ((arg != NULL) && PyArray_Check(arg)) {
        return PyArray_NewFlagsObject(arg);
    }
    else {
        return PyArray_NewFlagsObject(NULL);
    }
}