/* -*- C -*- */
/**
 * Author: Pierre Schnizer
 */   
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_multifit_nlin.h>
#include <pygsl/block_helpers.h>

PyObject * 
PyGSL_gsl_multifit_gradient(PyObject *self, PyObject *args)
{
  PyArrayObject *J_a = NULL, *f_a = NULL, *g_a = NULL;
  PyObject *J_o = NULL, *f_o = NULL;
  gsl_vector_view f;
  gsl_vector_view g;
  gsl_matrix_view J;

  int stride_recalc, dimension;

  if(!PyArg_ParseTuple(args, "OO:gsl_multifit_gradient", &J_o, &f_o)){
       return NULL;
  }

  J_a = PyGSL_PyArray_PREPARE_gsl_matrix_view(J_o, PyArray_DOUBLE, PyGSL_CONTIGUOUS | PyGSL_INPUT_ARRAY, -1, -1, 1, NULL);
  if(J_a == NULL) goto fail;

  /* Numpy calculates strides in bytes, gsl in basis type */
  f_a = PyGSL_PyArray_PREPARE_gsl_vector_view(f_o, PyArray_DOUBLE, PyGSL_CONTIGUOUS | PyGSL_INPUT_ARRAY, -1, 2, NULL);
  if(f_a == NULL) goto fail;

  stride_recalc = f_a->strides[0] / sizeof(double);
  assert(f_a->strides[0] % sizeof(double) == 0);

  if(J_a->dimensions[0] != f_a->dimensions[0]){
       PyErr_SetString(PyExc_ValueError, 
		       "The length of the vector and the matrix do not fit!\n");
       goto fail;
  }
  dimension = J_a->dimensions[1];
  g_a = (PyArrayObject *) PyGSL_New_Array(1, &dimension, PyArray_DOUBLE);
  if(g_a == NULL) goto fail;

  J = gsl_matrix_view_array((double *) J_a->data, J_a->dimensions[0], J_a->dimensions[1]);
  f = gsl_vector_view_array_with_stride((double *) f_a->data, stride_recalc, f_a->dimensions[0]);
  g = gsl_vector_view_array((double *) g_a->data, dimension);
  gsl_multifit_gradient(&J.matrix, &f.vector, &g.vector);
  
  Py_DECREF(J_a);
  Py_DECREF(f_a);

  return (PyObject * )g_a;

  fail :
    Py_XDECREF(J_a);
    Py_XDECREF(f_a);
    Py_XDECREF(g_a);
    return NULL;
}

PyObject * 
PyGSL_gsl_multifit_covar(PyObject *self, PyObject *args)
{
  PyArrayObject *J_a = NULL, *C_a = NULL;
  PyObject *J_o = NULL;
  gsl_matrix_view J, C;
  int dimensions[2];
  double epsrel;


  if(!PyArg_ParseTuple(args, "Od:gsl_multifit_covar", &J_o, &epsrel)){
    return NULL;
  }

  J_a = PyGSL_PyArray_PREPARE_gsl_matrix_view(J_o, PyArray_DOUBLE, PyGSL_CONTIGUOUS | PyGSL_INPUT_ARRAY, -1, -1, 1, NULL);
  if(J_a == NULL) goto fail;

  dimensions[0] = J_a->dimensions[1];
  dimensions[1] = J_a->dimensions[1];
  C_a = (PyArrayObject *) PyGSL_New_Array(2, dimensions, PyArray_DOUBLE);
  if(C_a == NULL) goto fail;

  J = gsl_matrix_view_array((double *) J_a->data, J_a->dimensions[0], J_a->dimensions[1]);
  C = gsl_matrix_view_array((double *) C_a->data, C_a->dimensions[0], C_a->dimensions[1]);
  
  gsl_multifit_covar(&J.matrix, epsrel, &C.matrix);
  
  Py_DECREF(J_a);

  return (PyObject * )C_a;

  fail :
    Py_XDECREF(J_a);
    Py_XDECREF(C_a);
    return NULL;
}